From efdb01a8ac2d5cef85c651177aa18ac69dd660db Mon Sep 17 00:00:00 2001
From: Anne Philipp <anne.philipp@univie.ac.at>
Date: Wed, 9 May 2018 14:15:00 +0200
Subject: [PATCH] whole bunch of modifications due to new structure of
ECMWFDATA, added basics of documentation, minor programming corrections
---
ECMWF_FPparameter.xlsx | Bin 0 -> 20494 bytes
python/CONTROL.temp | 14 +-
python/Control.py | 259 +++
python/Disagg.py | 156 ++
python/ECFlexpart.py | 1253 ++++++++++++++
python/ECMWF_ENV | 8 +-
python/FlexpartTools.py | 2370 ---------------------------
python/GribTools.py | 16 +-
python/MARSretrieval.py | 350 ++++
python/Tools.py | 461 ++++++
python/{UIOTools.py => UIOFiles.py} | 90 +-
python/compilejob.ksh | 14 +-
python/compilejob.temp | 8 +-
python/getMARSdata.py | 203 ++-
python/install.py | 169 +-
python/job.ksh | 38 +-
python/job.temp | 20 +-
python/job.temp.o | 10 +-
python/joboper.ksh | 34 +-
python/plot_retrieved.py | 3 +-
python/prepareFLEXPART.py | 55 +-
python/profiling.py | 69 +
python/submit.py | 69 +-
python/testsuite.py | 0
pythontest/TestTools.py | 37 +
pythontest/TestUIOFiles.py | 64 +
pythontest/__init__.py | 13 +
src/Makefile.local.ifort | 4 +-
28 files changed, 3112 insertions(+), 2675 deletions(-)
create mode 100644 ECMWF_FPparameter.xlsx
create mode 100644 python/Control.py
create mode 100644 python/Disagg.py
create mode 100644 python/ECFlexpart.py
delete mode 100644 python/FlexpartTools.py
create mode 100644 python/MARSretrieval.py
create mode 100644 python/Tools.py
rename python/{UIOTools.py => UIOFiles.py} (59%)
mode change 100644 => 100755 python/getMARSdata.py
mode change 100644 => 100755 python/install.py
mode change 100644 => 100755 python/plot_retrieved.py
mode change 100644 => 100755 python/prepareFLEXPART.py
create mode 100644 python/profiling.py
mode change 100644 => 100755 python/submit.py
mode change 100644 => 100755 python/testsuite.py
create mode 100644 pythontest/TestTools.py
create mode 100644 pythontest/TestUIOFiles.py
create mode 100644 pythontest/__init__.py
diff --git a/ECMWF_FPparameter.xlsx b/ECMWF_FPparameter.xlsx
new file mode 100644
index 0000000000000000000000000000000000000000..21bae48a2bef0dbb1dc5d02d0482921cb980a01e
GIT binary patch
literal 20494
zcmWIWW@Zs#U}NB5U|>*WSoJ+(NhKo#g9sY~gD?XJQ?zq_UP)?RNqk6UL27ZVUPW$>
z!Xg$XjRg!$45MH~hrpSTlX-^>1ono%vUk}RcEn-ns#P2^yJVtOUa-%oaMBUF{fS59
z-#+2ZD^^9cIArN}CRN`v65f9=Au{|V6O&ZB@kP(uN{q9ubKjT<&wsrA+G39hQ?9Nm
z5;NM!B+~rvt#AIGPYs;m(PvydpUhu!WKT;};wzhodzvvTL_em_@LK<A=B%bS+ivq4
zl}|em)Gftx@%>(hh0jD)*L%msY0s3<p7i|d&NeOKQ@>9=>55np6Y;lu$<wYwPZq@5
zt(UZNm3PS0VO+5;NS?>IO;BK?p-7Q_dIQ&a<-0q)b{cUVzErU%`8B)5thtBRnU%`h
z^d*?6IE%iBz4|6rPy6T7LphJ`=-BcV<z**zUy734eYLj5Ud~N2ey%om<{#Gw#-%0y
zJNtgyUvjK@?r`k+Q{_$rb+(2rs~Fb*m3F%jd}b-v0{0`^E2~Nw?Kk$ne9$4u_IzTn
z&Viq?H{9HU?|$N8Y2<n1lzji`ZNnm#xI1MTMq$nO4(lc_`@Wx%f#LsuW(L$8yEW86
z`3nOBgAWq}12+Q)Q+!csPO-ioh#bwSBRYe^J-=p>um52Kp0@WMf7(xczj}9qfRlH*
zqsMN=kh%uXJN#RF7pq9$-LEbq5F*vg>|VdO?$4>)cWruKolyMt$Yf2okchSNw(_-Y
z(|;9~vTLaoO>-+;^YDmZ+10OCzs+>_s?;dk+po2H$BKpCC+{n}UH<5@QQ*8@xD==G
zQkga}+nJ3~$9yC2CH+{(c-1NR?b@~t+szlentXnC{q<?fUh>*^So*xVG*k3MaZ>ut
zbu}k)6u<tI3JuP%%-oPPugFTNGCn2ez-O1)`-T73==#k$5u{`y^uy@9`*GC;U+bf`
zn0!1vxAKTN=i${SKAddomS1rXYbLqJ`B#I9k%8d~GXsMZ0|#S8jy@!xl;;;^C*|j7
zgDN&yo?+Bjz`zI(yb;RSGx>bhVFR9H@0<5FU#u;hl$SciYHyF@`snWrYJra~95I{r
zO+DX!vhyy7lX?rqn^peQ+Fw}y_u=yQN5t8(=7vsC)#TW+bam8ijUN_!?QEJAu5G(C
zck_w_hbHI3`t9|H{fhs%Ec>dms7sZ#?dJO>OD-j>*DRFFopVL+%dBKcsh6dd*LY9a
z_-;4)bG0F9nP%VArh+wQ3VzM~=l*B!-j(NFQ}{H!f7<$p6ZzujQUol!UR_Yoo)M*T
z{G?jv5ry@KpP09uoSv{r`f!T3!EMvT6E7OBTDQ+|+O;O@T3l0QlxSPKURxa(Pw&@B
ze;?2KCgabk|1Mecc=@gNId_vja$_rI4@hO|L^Clkd}m=`;Adc9sL0VrPCIKteDiM`
z@a*NkQ7=&9|5tFgyxW^azD8>fSc$~miYp3Cu@Y1Ao~L5Y+4lGQ#cy{iUoSmewoxOh
zrm9B8bhk~&B%9`qPOGP|soSY8*m17)*0L{&UlWV`I!k??`ed<c+BPkne)D(v=DlXS
zd#_nHF`QbqrDUC-&G}D4T50|zyj=f9ZiSqAcrC{^S~_8;E@O=9-V5LT_jYl*R<|74
z!|E*F`e4I_Ju4k{${lG);d^;I=t`Nx>L1Y$w9?9aZ6EGhC77hVq~JK;lu2Fc=Fv-e
zZ+rOs7C6FeG;93=?SsABb`_{|JpHN?|5W-%oWK{&^vzBl6>n12*Q~SOylSs&kkZW3
z(=iq&f~&fs_qN}cS+;V$wBpBn0oR58d*3V-DQ|lcu)%GYs90jd@#c>vJcqxjsFkY9
zN*t9rv$bf&L;tV;L)Y_F$7QYI+u6Rp*Py}o;;Ct8Kb)>OWcICZ=Oh!iWhL|0$m=ZI
zzb}aE8ute4<*&bbBz@aeSFdRFICq=evfIJ0-r5JBPs?St=elt_UTIr&QH7A!_1{w-
zS;)A5+gEJgd1tD|>z9!ou@`^aHQkuTH`npS;>m^T^AeAmy?D8zC-0RSL;jnvMY~Sb
zXZjqpc(7*c<HE~f(}gB;n4Iyw=h)MK|CP*=6nBHg9R(s6L#t&PLkyUfE&iOb|Bk_1
z556-NlVxVFjJ47#lwg--x7gu-E!8+?W1rH}`%NeITYc7F)VY7ovB@v~D9xGUe!0M?
zU+kfgU9{e`>mR@Wx3;?Bck;upEjr$5cTO4XJ8b@1t}8pOa7xwHiR#np65jF#-UOvG
z%<8_oDKBC&8v{e16sUm?PKPBKsky28VA2p&j=hOI-z{w}vQPY%Jx_7@pX{_8{f}(b
zIqA3Gs42c)@?^!<)c#GgCQ6>%ZKxP9DS2VQtZ8vvd%p(UViS6y`Zs%<Z)Wg`{TI~F
z-cVlqG;da(f^Ottqs2eJe0jOo{QhJUas9i!D~nr#{UR=vtrpi>*F0HJdUDFGvV1=F
zrORFyG)^*Hq;oL0ym4oB=I_7X^%m(0K9sP&VE-_B#alP0Guuw3+BKye(cd8H_uyvB
zQbVDF*?y-Vi^w+LoiQzYs#idtpl6J_qv4FbW-})UmbG14vU{(|z2yb7IrOC6Ozuu!
zw(+8K;_>#-&IjCbm#$bkZFzghBDu`wNO|k)fFrCn8M{TlKA*N-f7->cd!1dpwOi|$
zzfS*t+Pa57JL^nI)ds06pKJ7Ata+o&Ao_Yo-D3GK%~Nez9Tp3oF?qb`_Ltpj^Uj~&
zD5bG`%Z+!YRbN*%UVn4_x<z@-PvZ|<Q*M60_MzF4BhKUl*OZ5~-TqQrFMxy}Hk;Ng
zzVxz^mwCN|pz6toQoA0pe=V@QEt+!ZMW5)sPkS8~Pm|UUn83UG%_L`^2qm+pQ9bWA
zlrS)GF7cRfM6g?O-U7jt!VV2fkIv2B%LG;K@u>DV7cUKy40-h-{+~_hQk$l22A8V5
z;|{c}H}$mRlV6&mcI#bx!quw^`txi5f4V9B{ORW7>iqNO+5C9<<l*7T?DlEfditz2
zX1<b*ulap*eti9#KOdjW7VqD;;qlJC|6kkv>wf+$wiZx}u6X!!dVhWOuOC0d=l}b6
zSDAay<VDP~b}EL?&V|;zxszac#<%idN+18F`fJG&ib-D_1aEAc@vi#>LtLyv*5T5P
zFZV=ia`f<4RMs7TuuFbvv7?T1d}k8#j^h$<H(hzvpk1VM_VD3Kqpke&&XwNbZRyKN
zZr3n$x;_8ALuBXK6p6&<XQf?2jms<!9-4UmhT^elH}9$)`S4m!a%s=z&VoGQ<=q`u
zRh^3xRDv#MwjU9g-rB3sa_hs{d%DX{ugJY-Ta%*nu+`y;Q(ec(hsRE5=mq2$Te=ut
zmEG6vlA^qGzlXcky7R0qRb86Jzmm`FTzh`2$g|a~f4MrkZQoRv2+lqcGQoNF+1Mj9
z`rd3jY~b|q%&D0w6Z`%)$UAi@J+~4z`ypq1U!an|(Bb5<g8NE0IKJO>DiAIST^^#+
zQDfw<yKzs?%%668FIW<n8Y%p4n{`~pGWeE%=dS;LS2nLIiMnE~xp7~V_3|Y(!td-g
z9@v`ruS>2-zs7yT`fuIx`|tk8?Kj_h`uX}>tpUw#;;T!ymi3$K-BoY7RJ(CO#q`xz
zt=m7D$T=*&6}DpU(=2WC?gvKqH6(Uqc7@2Rr#@0vTJ)FWY(i&zsPsk6^S=H&jgw`<
z{;i74?fj^DJL|te|Mwdfw>O5ZE-IhZ$?m==WBrfSWjAs=cw|kwVsD(<bL+a_)x1hi
zhi`w?Rkj@T5h~^GJhn-2(T8%CD`(ENx}H9?Nb{_|^VXh6r$ilcMMWd*EjPT%+P!c3
zT$TgNxB_ja6`tZ_opy2RLCfnN`9*41jTyv*PxELf-JKL;?&p6_&wFcl#}-E?Z}Z75
z)7k{39@|)dT|D`7Pw8rom5%<ql75<rm`I&97tGdFTE3-GZCB8RBbIr)9q+bZ-Y7V0
z>s`(_e;yjmb9`5@^POu6*YR}aiD~+OUTd7zdAg{*Vv1!<r{6L*LqUsDjn3Qa=KRx=
z&RAsPQn>n>gTTG0BHam7xA)4h+~qNv=kvkiMdJ0SCCfKQ2mJ8SUgID7{W|;Uc{zL!
z@AV%$evWbOa&v(tGw!eOV(;T@XI_x??5ts<@Y$&oM1PjAGI+XawY=sknNzm}=elaF
z*Vr{9y+C!5w$-b=r_68l<{Mj-+-+qKj!gJpdewG!i1+!-MUhv0oI_KtZs~qGTl<sL
z>O~z=Ym8oo-DZo+h|OMEI$^~%twOh}wXIDH0=3zUS9l(KZKPvroo~8w?&Bth&n_Dd
z=C77Y{CebdSlx?qtL)Uqh;`bbk5>E3|GLguR=|E_G5hla#acVKRBNg|1a*FGn*7Aa
z?1AV`?N6pkF_xdFY`I%4y8m&jsl86<th)VAuI}|(`Lg`g(!TfCSJ&>3dRdzNec^@<
zjdj0nrcY;_<9^Zg-$&Q2>uVc2r@QBgyXoi|e{Ed4reig~!PQ&O_zn0NFYXk-!nU?*
zwe*yzujdv#`MT?DkT%b*rEvmgRn2iP%%husuPa`1Aa1!_mE-<Mu2qQ<pPIUs+|QXR
z6}R@<N2$4K_g~*RKi$yavyWX-d_%S9`nUTTu(hS+3vK3bb1^Vf8Zj_PAlg#J8L6oy
z#rj~%2vm`+iH<FPv_q(Nf9P*{jeFu7elDAobJbyn07v1q?~iXP`aP055OaIxnKSF}
zPrJG)TKA|%v^9USlApZi+Q(tLb@l%JD7>+K-|@(KOOBu3cKzkn=#Pcd?NpCXKmD!#
z-`+p9m-3#vmeh28dug~Y=YP-VZ}Q>$|Gs`?dv@Q`^@TB88XnGw-xlBf$mDg|=eyZ8
z9iRNFy81tyEuJO2i0R*%xLrRryB9CN^kT#NUwiMZsP34Q9Bp*QNWAaMPs^J}yGy>8
zb@aEYoP1Opa^T13>rbX=Z##Uh^_@WK{d<1<-(Nb+sUbJHHt3Iz=k2wJ^R`u4E{puK
zqrLL~2eylTh12ULzO<)IUt0Mqrlet=Hhb2&%3E{G8ehJ8q}eFRv3>idA3rz!JvzI!
zZvWlAyEJ6K+*hePfAV>lpl#bnyL+M^h1Ue{{;KWE|4>gTdG>)tq4jyKhre$MTyrsb
zR-JH)w!v%1?UkNNyL!u9KiSSZUhPo-k!{1d|Hmit95L8+)#63E>9W+F>m09D^?a(B
z_wN668`1R!1`?-&|H$lp+iidO`)rHKO5Fve9P@<v=Ggv?Ir5V)jF0W|$M!|@V|TF!
z{n%gX(ZX`@^Tqy8zw|8GW%DmcGw#@0eY*S0f?ZC@na>tW##}t}y1eai+<oD;%8$Vo
ztL}R+-k)ZBSX5JO?#Ao0za2Ma?pc>8@Zj(y#&stQeyK2j2{Ael>bp2R*rdeP$$7`3
z#Lp?mJKD7$9iQ81FXO*r|EzTGi3?@4W3;)p2EF)@)}$$PK$P>3|1R$nzN=X6LYimr
zv^QToQh9ZcSJ8)sQQK!NIrw_sq`5yHJSq2I=~J_zV&6(*RsF*CyGoZmQOP^tbve&C
zxam{JBCCHXFL=-2?V5LK+K&L=jk_k#N?W1w^y;iv)@7@XJh_^+sLQPS?WqI7ZaZB=
zT*Ad8f8A^H<p`NDUuomgpT-B%`eOt(u~aBuXqd9lXA8@grQCBw4KrMn#TIU~e;3lY
zmTUSoyWawPPV4RGns?f)gW*ka#^1W9m500Em3LgWK5|E+hhsxhi1xjL)6&=WBq~dv
zxbiN=TRzQWm)y(PeC-LZe&5UB<axz@dWqY<bAQz36FsC%{15N_+Nmm%QOlXx?WR-a
znVze|`fHP|uDRAmj`-AO4vVe-XZ+67>M6Rv_H5VM@3%JG`f%>O0F&MC4WFx}+e$?X
z9<j4O+AuTX^Tr#wEBjC1W_h%GPwVV_Apvc!I?W7|=cn1u{(Q{2J6Fr+X-)UlBbV1+
zYrEQ3Ue_)l5cag@-p{?Phj(*Xd~7~?;`2=gXYCc?XPt`ekNmDb^tZ{acz@1?yY{K?
z)*DRyz47S!cMDCo^{kr4`~DT5;+328j{n+aK7;w(=Y+=}BLlcJ`=5E-_-iMq%@Xx_
zdX7M+YwhlDGxz66^z(kW)|m0uY0E?BJ;e$8$`7#2`El@|M$WMeZN+(sxq?dz)8o0c
zb6r!KJ?^X$ooSGMKf++1`Gx6k;kq_QZB!_K&%|Z%wfXRZ!)ZO+)FwpKPT+|QTY9Q?
zABVBW#Roxa9c7Q7@=4cIw0k3=;+a@}j;S)?fN00sOIukFJ$e=7bb_tltaQgN!>GL8
zX{OtHUQamI`JU&*-ity7i`kilE{Ww@9+^C4qnw-#!{KWecPK1XOz}5aKGPuoyl4m8
zX3o#w7qLy%%$>jX_14~&r)u7<n%WXZ(Z|$2UI=-ec!p^ucSh?1tE>A<4WoBAm~mZM
zZBnp&)(SUUvk$sWSugMGE)iE)5%#XRJF<dpiFTL0`$?N~j)6ZImTGqCeu;S^sFZnK
z<Y>@KSKdioNt^a;n`FSRFWTYdRaWN~#J50nQn07xjI|fn6f#J&E}G=@Db3SUVCs~^
zGnn4=d7V0_V`D5D8}z91{GAiBDy`>d&suV0n}zsFvFRdH{(te<xFRqpILCX(wv8L+
zP2+M_RkgNP)?~lHb#mnrnKRci%G`oZxTv=Jo}DyFK})mUkkvv>OOt)`hch81o*#o6
zi>DltYE;uyo3`_WsIRA?rB2qADI%ve*XQcAE!m*1rFpXi>_EN-uTuwedehQXT75!-
zZd}pP<XkkVNFul)dC{aI1Mc%auG2;Ig*ZKX7xdJdO`YT_q&s!ir#IVb0uP(ntv_zI
z*0on?lGNt)>({Etz0O)05E2}5JkxKP&H}~L-@-inOrIXu>t^g5Jg0iGLr5@R@~uU(
zSEn>bHl&3F^96I}1}pxYyG}DY@QkMI^riB#!8zw=bzf&rn_;wRC0A(MN=5IJ>vDDW
zwCoO1Sf#nf&ouiK;})l)RhsOZO=r#R3_jsv`qW{cD_DTx?22VFXAT($Uo=f!{e`!U
zuh47So)X`UzS(;Bq-V}bniY9PW1(jK0h99F<e3tZR%NcHt)V;9d?Ovw=1fvprOE!z
zl<$wv94EC^n(Uu+Ck6Xj&S1@OGkuyMc5B(E#5s>#{WpHM%=DV}xx`m8m4zkK^y0ig
z_Z0y_!N;DZCmNczf*fdJmVJuJ$W+xaB$)3`_D8!=k^L_i&g|s}@t;JiET2_mz<r@7
z>(oJ=Exu2keLGh@x?lBb@>7r4U?1J$kE?F>1n-H7SiEFe#IM@i<d2dkgBn*X+mLi?
zjVAYrD^pkkeJu;t7y6`|wEBbu^I<d3N#nZ6G0o}UdrBRnrYvf$dr|!0%8D>W5B`{(
z(=i%a`)iiln|;#e?v&gv;9}b5dOxi<(qY*#m5!`4jn%!Mtri=mERrggxGZ*Z<uREv
zOeTw^if4Fb%s%BVC7iF|)7v=7UVX!)HtwFMp2e(}SdV@B<Jzk=pR00Hs<I)=G6$E0
z8&_u7-@oC|c0j;|<=LUVWe@LZ2nf#TzANVFlAxEy*wk?#W|`1D#+T1K|K+vDvS=%(
zZ8XVT=df!@L~lai;wR3E{@d2*$kq2R3kV6B_Q>%1Y2I6=TKTOCT3YPywsL=WYh1BH
z=EU|4ch{#0VhbciLsOUa`LujJa8xtyEyJv}b?+w4?Y^$lmQe9Grg25WvOP6%%}e4c
z)&|%5$8Y+a8m}cVm0PUsq2kOrc}Z5g9S<Lyo6)MUN?iQYR_>o}jVlslvTj9hS=>~(
za=~Ir(H?KsyJl1Qw_C<J=G*CWrr+__oX9*mYO=^0n;O0TRJliQmf9?vx4rD;Yzr&9
zMwXRXQ;$D?|3^~BnC0N+Z|w4YSEFz3|7C3PA>+}0zO(|S^Evxv&phav%W6O8?v&uf
zzxVjk3v|vj$)C+Twb=V^RaCh6@)xT6x;IUKGGqI<kH$WNr;ff2Q|muJeSY_skZpd)
zn&+9z@Gtu$Z*#zY8S9*`{|`=w|2XQ(<IZ0#`{8%-EU{Po60=^I_H=I#el0Thk&aK@
zLCMM!O_y{Ug?~Nu*V4Pa+(dhmhUTrc;Tiv@9uoh3RI_n;RlU`WzGMGHIv%pCxH8$h
z8SFcecJgwnb+m2eQ-7{i^1U&CS+O-Bt#ggM4>B+?I59FXC}T7rAtSrsW+Zq#1U#<W
z(toh)kb#KncjkZ271fikwemOzn=9=Iea@Wp>lWMX=@Y(-{H@+BlX~Fk(}MF)&1V0$
zSt*^yH&N%xE~eHe$>JFsK2+b*-S4>3;?*4yPF<S`iI?^E&0_ZR^Xfc4AwraUV{4$w
zUx~wYM{aKH<6Jqx?%7Sx&Ps!mKf9$btXY+n`zfA%i@x91_c~8%!<(&agChUE>8q8u
z$y>btoX4m0EuNMW66_XvSe|fw`>7$x@pbvPdD*pmEX#7woc{cc6I%!fAJGv$L`(>P
z8`7{4$U1Bw;PSrpPxFmx73=nZF4cH1#p`N+8C>ITFZ$+{v0e4&?-w^(ocazY-!q7<
zj}-sf-I}b#y*AH5vDL(7mQKxFyKU<iMA-1wdA5kHJ8@z6S?}+i-P5BUS*C?ub;?N;
zJDGp1U`FBBHNkGHQhZXs1|`q*`S$R(e@2<_X7i`(neWUEmSy99YOrO2-}RX3^LJmi
zIM(o9^UoSk5S&z<lwi5Y!}dhz+tV#cj<4On%_#pR!`but%-!O0W^6%VvXFBUQ9)n~
z>Pw&oK_fT_)LVm^p2j#STsQm65E^${2^0iNfBt?IG25vh6a;TyR?oX4F4otn^QM%U
z>*R8~pba1DcE5>p*l6+U?i3DPn+XrE<eI<L_Vdf?JU$~b6r2j=dlT&+y@^=f6LNCG
zsW(fG3Hg-$d&e}(`LfOIPq!KG9$L9>Z)?%p?+fz!W2XNvU!K?@cm8DkIgd}rTR=fj
z={&J&!s=U7L8-tr`nmdFK~A6RvyOg#$Br!sK0kPIrk0g~;g<x?juhDQb0TBc&$bY#
zjnDkg{AjPD#EI!nTBhF@b??$%yKVPdV@1Cvy(2mW`+vPI&rzu0Vabs95|Zyb-_ve<
zto-*++s>`wvz`j>-nGke_G_uVZEqxE&#(LU-u~Y5sq3D)7TB^G?hXq7zhLLT%J0AP
zZ<|lEeJ<|B%`35P=j!+_rMnVta&Z|#PmHHs((aiPF^z?_=IlAQd1)pVU9%_MQ26}w
zkVn&tRnvK9=kTyw&s{XZZo>!Fh1VCUoXn6BSzy24B{Ocb>O+ZLDktW}U-DP<4q9O0
zQW#Rh5ftlvsJDLKCLN0-0^fIMvL~L^-M8h%WQ{Ws=lqWHyyywauzt36qgg}dc3U?A
zM*rQ@E>2!uuQ+@Cw8hf9*Lqy|-m&lLrS#tCzYp9nmjBSd&W^on&M}iG6SEm(;$F@>
zwsH0U#915C{&R-?lho*yUjAq2v+qZw{!YzoS6E&x_P?X=C~FsE)<X~fgWfTJ&fc4s
zQL<0u!)oTp2PbdH*vn<e9nxx;At`vtS&Uip&w5osef1MbAN1DctE$hOI5pXdyCM6b
zcie_0^LkG2x$Ch?`eD@FjPpL#dU@6t_IyxU_e8HXH1NRIMGbGVl=HW(p8YWDe$DLv
zeK|86j&NU(pSpL2#)GKyd*(^~+VrR)`)r|k=!XvOn+1PEZ@%tgT4%~SbFIwrxHZ$7
zGtEzJ5;XlKAaXEo_C|BwHw@G4{h9OhQVw_T+Wl@pp2ZJt_uW=1d|Arv%O97f3uR~T
zaXQ&9`fMX}@>;`ZE14gg*sx{n=4^R1OF`g5uZi^eUC%?-<vAYp{3Lqeu`OFZ&*96q
zZ1TKTXP;=K1z0GiiJW)poEf||scees?&A+ac><Dx(w?7cx|jUH+Vo2c-{l_u+Vd9$
zWOuCCvdD{T4_C|dSkK1z#$StEPBVYCcT{%Y9+uGg`FW^<&E8JOI`v>-j!PWJOhDo$
z>qL4*IkHu{j`;-$b6h&{-FETb^z#?jZVl@>`kLuxm-A2275V~81n-nJx%Gv8&d@yH
z72f#aOBe4?W`!ja-z%wgh<uplwe-dwKas#!*W{fN=hj>{sXKU7ZuQb&O^fwz4@*|{
zo4cEeJ^W&0@@z@o4)G7u9G2eNBh|{PrS`OZ@f3r7VJnK%JBqDW`oBB7&M!^V^+oBT
zjz`>ey|1_hUeCF5`&(N0`wdgJRG4yoW!xRQdT#BCPPW%`u66vFw)gLo)pJ+4J@nc7
zsk3cv8w>ZGE3G!0S6ANXs{JLSYOChU_VG$lG2`N0`_oMrbqlq8Tl=>OFE6ueU7UaR
zNt@A<?0lwO7q71SyUKe(Zu8P@);^C+ejeeRy!*rk0|$>oQAh9HOt~VLUS>97-Vwo_
zjXbX&XL!!vHB~U>LD!D1scJ8-YfOn-DHr&!<Bmn!SH<<)ujKY9WWNm7Ug!3(H&>YT
zn)>%%#acPd`*l-TR!$9F{Oc7z_k_Fk5eFnbfA>^h^<202#5Oyfgr3#En?Cy~KIA!3
zXRDLY6Z`uQOKY|8x+n9_>P-Aw;nARfLHE@%wt}4CKi?+a`5ZmJf7{<}&r^GR3`6=C
z|Bhc%>v-*!?$<X*PKX@;Vi7ef^6nqbP1nxY)=u5;bI)vLb8Aby?b2@zVi#sTikEZu
z{<q5Fx3Jv1lk-haPPB;1y0}XGQR2$rM+(yXQ8{hG?fgz5-#1^_u;t<5V6%0`X<j)?
zZB5V5|CaUlc!5%W#!DUk*3-<}O)M_kH@}gIE!JJ8e&Tx#8-H@xyq{g?S-p?0U%k}U
zeD2b{djA_Y%S14;@3d_exA?TYz{Y&O<2An?(XZkb@1Jh2TYS`e)lx^V&G{GoEi&7$
zCPuOxKlRAx{?*`VqLBikyUH5x*BeQ6$X_~fHFwjt>)oo6%!Yqh7k3=$Ys?XrjkW#U
z|J3rw+x23*|6`qrKDNI5(+eI3hWQ>ii)?VCWlhZOqDNL@wewf~WbgVeyy2|sl1Xz`
z6})7TubDaTw6X7;7ApZ`hhN{mDn~|d+p@iH&+J)s866STFE>q3$xN*I{rmghb)UCd
z%-i<#cdYr2SLN16e=n=Kw{QQw(*Ix2AI-mh_uR`pb@g%YPkn#1?)PGg()#`Xf8WpB
zKlQl1f4*zsJC3&Rmt()bt+BiFrE2-v=#BC7*PU|ysApHH+5h`%<e$o#w{!W9C;g3K
zeOtbF!+rU$-(J714mZ3bV=43W^y#wTGpEiT-&bL;Jt=u(^nBOf{`t2Ldv1PfF0<?J
zuOH8kKiuxS)1ISC|MuN)e-<mfE&Dcm_VVlBKkfLpN%-$SIi73MXJh|5zW9Fo>_(fZ
z`p>l;uLtukt)2Hp&GukP{qH@Zwi6EDvlP=WzuW#nt$*F_y}xf~HT>JFyKipK_qp>^
zcmCfKd3sw+x60LewW5adRRPbQct`%zY@R>w{jS=;vny|iN?TgZ-CFT8%BOp-!+NHJ
zg}?GcH`nN|dVBTjtoNnTVxe-Ilb@gd72<RHZvD3N&&?Gj96!hGSz}-Jynf1kx$d{G
zZ%0nyJL++gad+iELv26liJBK@tS-=aeRGPnqt?IoZ?*>&B`vvrfA#hJF5_KXY=0**
zl$*&fFWz}Q)NunBx1!G;rwD@`QZ0u+mR=B&bX-z#b>iRE4|m9KW#v-{*yHq|<8si&
z{#DJ7La(c@`c{(5!uV*boM71bV*(m+iW_F_v0svQ#GrbW$c|OLk3^N&!ykS%kl3~~
zk*V>KD3iLIm|wcv-5@6~#z$w@ms+eis(kN--sU37jK3$$OP9)svP9W`^1YjWErB)4
z@zYuTZy~Ea{MPPTz9(wyye78*=Dq^vwdz0FUdhDwHvdzRvQ|&`_kEisJ9F*%*S~Kk
z|6Hp)t<><?^21SaDu<*!)vkAB6kXkLg}c4@ZpRU$Z4+(gt4!n&y(z|{9{A_09Y@>t
zz3Y=s-ePzm=#aQI?!*T5piLpFyeefz{W;I<vnIqGQ~BIAWyvZRpC!v$7hiaG(!}lZ
z6E{cyo=DpQC-rLwc8Z5)oN-$8e%It(_veQ+PT6r$Vz;ip9*@+;rH5Ps&&n9j2<Lgm
zVs4o#@VqW5o1xWj<)<T7kL5iiQ+bzYO)%akcjBX|ln47Ws||+Ut%mHA4(>@ZP~UNK
zrLqP`V0DrDig@!|QeGRYVhpAjos#kjP2hjZ!mFAv>*dZ6)n^ARE1D8=gx7KeHr<}Q
ztHQ}7v3T-vDUppc{vKkPnvrb~(rD_qe)hxl(`Kl%t4S!tCO(?sxF<<KeTD_ZHGiJ4
z9o=K}<M~pT=*P>sj{JBz^-Jdzox%W~35Bv-lvFwSyf{_Nz1P@v?!B<ye+s*#>bicX
z%&$VtQw&%;B;p>w62BQVgYm>32F0ysWCgYu=!hKIF=Yi)@e<FILhMHw4Bci;QL)fa
zop48IquVa4BNKXEZ-4q1(Av3CY|mGRMXDKG64_acLKCkgtm-|%vDU6woPF08J%;3-
zEqy}aJJrs1{1S9{<q-O|&?3=0DEbU%aJy?^+Xquk#XFkq*2=agqYve<cC1yHnZsNd
z`B_he|I-77^w}#L82rr{B83iYUVBi3bAj){J4YQ=7Dy#~IZUlM{#$T`FQ>twX#&@-
zh`l!ZwQ2XG2BAZ9iw+rGn!HE#%(RBrB?l)k$S~?~E=VoVD%LMa@??<ZO5=6&WLR4M
z?vaC0B*U>O3|HR7vx=OU!l0bvnbln?{Z}ZVK)uEKX+b(eR@J2KEjNV(q(6IyF-o?D
zGj2KI|I|;^n_=y-&1{QR8dQwR3k5h=RQYf^oUP;3n4rn%c0(vCWgSz`hN<s5{yQ~P
zPbxZO#4uCngZas64K^Z=);d^BoTn@iU)S&W;DN+nh3T<6+%-kBm`+4VGEKQ0$-JsW
zaFOq^JDyDr+zc}hH85|8m*7g_Y8Kq4vi0b$iT}1X9-J%iZpz`gj3-3e_M1gZh(|K;
zUrJY2Qm|&Ru%FGwZD)P1!Bb<Y&_UCiN;^C**d;TnEaOl(qtC=@<a&T3deW_=yPW$u
zB<^rN^s~Gp+2AAiBYU%J!Nkgs3bMgUOeM2sF`cl}T;O|P4Wp*ugvCy84NPYLo6ylP
zLvz+OeJ0VgE{3ent-6i3_qZ4w3Y(**+G72*C6&R*^+En5U8bCCA!a4NGXrl*UW-(8
z`8nGwT>8C`)284xj$SM&vmYoaDk-Q<T6Y^}PIrTa_ncc{G675r+!8d^b(kk;GS&O%
z2qql&>}n{C3drl%W!f^yzxh^8mRg~UPqUPYYeTC^dEtfD28l$goPgO;{Q)UyGONC&
z&-sx$rTVJK->vQw8BUzZzsvCUdg_v?mE7Cs`eyPzjFnfO{h(oHy8Rx@*=$=6rlw4u
z^FKRas&7$ldBKI8B<?qp*WEg^jA`agmb`j1!*ZXU;l2~~XCGOy*hf-iX7T0wUfQ!4
z?on^NRr%&UXRmFIYinPb;Jtl}C(J4?JgvRfu*$W%=h^2wu1-Nb3dz5h9p2Mz*3h%0
z_|UW^+hz!FLrIQeZKt1hN3FlJ@bJ^_4M&+3qFE|_Uu#&xp_H>?x5b9j>=MCDGSm4F
z`EVA**abzex&7kq!D0tcOeIRb^t>mO=6&GlCqtD71wWOxyr0#^H(mdboImrbYmJJR
zF5dBM@|9zlsr2ASjX+vg!{g5rZzbLJ*w0~MCl8Kku>S4?Pd^*FJoxZa=})Wwiw{59
zUUV5|J6wHgcPMUscZ&|okr$T(|79~8{<28)W{_NN;di5p*^8t66z6G%lpLX`Lpi)2
z@1I}a|CxJgYD(?{_9syc`zF17&hvpyZpN1R&7RMG6@N~f$Y^(z@k9*ApLt=W`UjL-
zW^yw^vX$~*ZpLP7q?|6e!)Nce+WjWY8|NDzsN-u`!+WqHha>Cqx?NHauBN<GQsi!Y
zWMF#W<Clh8R~U|8ZZ;@ZOlVnK;$GwA33ii*|EFmN(ad`O{R<h}mzyhCGquZ_A4#{b
zT;KMo=E`=v7fFWscKLmw&P%qjHuwl$nKsQ%#h&SekLB@$%6B&A1r=|IF?F4HSU%B9
zGuSO^=JtZyFZZ08%Psrv4#VzZhE2=9Oz_FPJ$Lc!ns#gPiM(=!&r9wei<Ij+Q&P0Y
zg2Oo7Kso2C$${w!civ6eX0818T=b5&EW78JAA5H8;U3-%pT(@^Y**YCFFnEP<l>(;
z<&zXNbs195g(qd%GM(@-I9>&E!X(RhNgqzA{m}~6m3oVNC)%`MG*ht7Z}&7i62J6R
z@9dp7CMjDS&tw)duH+GZD|iTGzj;vf8taGmOAU7jADUaV!x&uH2{oSi{{Q6s)hj00
zn;-cn)9@x$dS&{xh9_%|m8YJMmHAQkhqZwx_JzUwi98N%de9O$x)3>my?$g=?f0*T
z`TentYmHbZENv9Dn`W-CRr>h-byl7icIGiE1hdFA^BwXLEa0&Uie4lAGFHZ7K}*5%
zIU9PusN{2Qc73pq_r1_(yP_Fo&NcC34LOMmnE$R}4Be-@@Bgurjs33}j7krFe4}vY
z3hVLB%Ag{xF1EnN8sUj+M^op`{jrigQpWAa3O)sI=HyKlSB}MPwMfuTd#5D%C|AD0
zN3h`Cxn_f4meLt^cXc_-*UWz6Z)sp!eqFsq^Qi$S;|r_ROXU1h(%p5YZuOxHU+Y@@
z5`r0SzEla`_3cH7)yMCJdP1}OT$7efPTtSBedal<HJ<9v`bsDNdC8(?Ja-1n5VuQ}
zGbf5!_G}R}3uH)HF1&jsug6m1`@)Y^92c-QoGB@qW5Qv4IkMgK<q^GDkR#L$X1+Y)
zquH*hy!Om&gGuL3+iX1Xu=8R3W*^3s<-*lTydHNSwb&Y}b4KQztj)7JvtL&I`%cqD
zUn|4%GY{i3C(QIaGo@&T5rgFYx`?$~#G|LHU&||W<MYpH(G_O6`SOUly@|q9zduWt
zPi|Q<xgyEx%zVG1<59xbmQ21g{}{x8M-7EDj2LFVT;pqLWUXI(;7G9Ji!hGa8Vs8*
zcj}&FTOwny>4sI#7Q@J#lV{UsU7PPGnx?nyh($Eh<ix_F9pO?b%Y?TZJvB&WFuI(%
zE$FP@*CirfTbgSNYu3(NS@bIB6#vgT!auh(?FzBI{WQ2jd7h)flZI2lAJ2$~SIXS$
zQ0MmeY4YX!_V%yY8>8gox4&P0Ze7Jnxs^&=%0kY^yo~Ss{``}I{1&5qa{Rvz7TC;K
zQJ=K*SJ1qpi$nfBI&A(iVxE)DxdrvA-=^rNR$i;`P&*N2FaEsW`S}Jrk>lrY&wrGu
zbZC)QfAq(Q>92NLm*1UR;cy{S*-rd?nSlHsoulfzZzgKzt>sZ~KA9oo$FjNN%!NvA
zxkP6>i+edoFK6}npLe|a^4CSD@Bgum>FKX;o|MJHz~Cdrz`zR`*ex!p%t<W<4dI=M
zyq~|=LZHt6<v;c#)=%9VoJ2O=GU?rNY2KEg>PvgA5`?mQWg9R1ek=L^yO&ob<+hr`
zB;Vf`5{!HH|4g$<Zo9bL*KMo#9PyUT9^2pR#2(gp#+lgd@&9*a8q2OQjw3=YmWL<q
zJKFj$^Urhn#ZEbzbx}z=CIJV#JmR@FO3g|W)Guu7vfg)y=iL&U^Aj2p_Zf7D_5}u=
zXgA8Tn)Y<Iq>_)mDbMB|y_0yJ-oLJ`{Ol$dr;+-Uj<uJTW=Gom*_2@8Ud$C>zODOF
z)OM*oLGveCJZJSkD!anCa@Hh?{ce)87IMxyws?whV1{A5vEWO+71OkN_8b3eKKn#h
zVZ$V;{;v7AxEC2DzkXIBQhch|m%*69=d8h$he3}R1C7piig*XQz1TLV>`9B#=1HeH
zW~izAb}n`?mOH!jite?d_1V|YEnB^IafVKVyVJ{atM`>C`1H67#4<h4VPjsqtM1;q
z1GPM6PA}P1%bltwE?)b#;;p=fg=(kE#)=HK-u2l{Cw6Rpw@2U~o8t@X#DwU)1yjyz
z7akJYbI^0i?=FkYiETY|%qHA-EZFZ}z|X*Ro4-k+gnirVRS|pbnhu_|<j_CSSj2YZ
z{+#EWDi5sITRiw(lx4HXLUhCa;!P2!;tz&@s6W@*kYg&R%A#xG9I~wCSLd_I6Lu)N
zyYGr{>RVfUVf{zN=_|rARM&6N&{%tfv8Z>))SApYx3`75F{K(V`{F6}-GA3(-(~lF
zS8XgcJAFj_21CH9ltpXw6+f*$J!x@~!p`h0?NakAd?{`()`j0*vYYi>MU?Yu(Ytla
zQVtsFl<(gn^vc!p$}}5w`-5(g8A&_(PygEBdA{k-u|*D-zpnW)aq(5-dm68o*!LYQ
zKd>^ty6AWXXS`U7=k=$Oe|QC+boxx#(Hwp7Lf#&kj~VIA`YTuO-oI<Is<E@EOG37>
z>>nu?h0wXz?{WWUi&QU>Wm#eqHs`_}2hQo$Qy#ps+Bqe~H-zz}_?xKdXTL`a%ss^*
z-#_j9x6|70vd@-mzH`mev4&%UdB--nsjpA4y`K1Todj!b_tT{cv5N}giznYP+g`il
zOjX#zJpr7r=4`jR{z{m=qnJUafT8)Y(Vaz0_LqBQ@8#@!GJkXFcZ<qW&X3X$cE1gD
z+86n{R$Fe1m0p>it>Cu3@l&jBva!YfKDt)_`rCVqRr^Jq7yqycHw#*DJ#2bbmV@--
zi`teiro5?=o0a3@{%i89kIy}oSDwAz;&l3BuN`NGui5W&-Kw)^JDJ@L^-tT*_rR~9
z>spA-`#iT_|9<JunNZ61COqi=32t#O*Cp@Q??3fy`rFf{3*0|!ZO9Y-ll6Lz%dh?Y
zlYbZdvMBC3FpH(l`oKKFs>v%K?&P2RvP3MZ@9T&7$I=HR;~B0@=gMfl@ovL{4|5gP
z_pCi&|Mq#q*}n@~*5BICu)EW?q~`VJ^yp(b2{+qbFFwHWr!oDE&Fz`L-aWZEyOQ@!
zoc)_^skh30&U9`1-qw^dLn!R`s&4_sVslsyE4$6us($8g`ht)D8L*AAU)^%|ZxAa3
zgN6tLg9x-x%1A6qO$jb3%FIg#mrY9|Zs$EV6R6d{P_MuulF|6})F$)H4IB1)Y`&Lw
z>rRw<rOpA5<ZcevEKe`LpeOgOC%rlH)T5}`uFB-_=`)tb?el7?PN--fj97Tpq~jhx
zGn>}6&VbdK1;u|;3q%ZF@9^4jg=^VC<r@z!zxn-M?V45F##4UHS9bnMmbhin+H?E6
z6?<#^v|A~69{O+Q%FSy~H)AcB$>V<bVd9L!g-l|vopxpx82l30;MAGVI>%(%ivpQ5
zYGG&8^2$#4J}4J@qf)ZjGeUAUzr@3u+Z|Ct_n&{+od5H#xy#|ESzh<282-L~zxl;X
znWoeSJO=EMd~LfAIZXFvdL6_4^jFQ|hiQAIw~1*TT$5@U<~(urK8t(r3T^B^pWpnW
z=+W#x{TlACl^<f~dGr}LFMe(ns61)K`lP9wB3>=pC35Hj|3A@^b+Y$D3SQQGKfYnk
z%RVvRs@+8-^1hs8tb5SqnQ`BBgKG5>S^KP0f_|sII{mKL%6T2H*nyjt+rKaI$k^T*
z`DurEbfM=Y(djGq&)70E-zn(wpTeU$@3a@PoSbk%%-`X`oXAa2Mb(e*7hlEy{M)>o
z{Q)y>uoSuX<#VnzUN$}FtgHmveXhWwn@$s3#frZ3%&r!^6QSOj8))t{`9Y)CMG>j;
zD{kv~&R=$)bk=w2uA*&gms~um6*?olzA^jljyKMCk2#)J3O;+*dfJVT399Ey9cPq9
z?@YKXWVkuBYnF}IzmA(rPI{U|COE8{FU-xU#8EP(Qg2#zy{h<nPvdQ`Crx`(cXZP(
z#;seHek}f*>$B@pXPBx*qsrk+N0r2RS93P*h&*`USjXnN<44&z#HNT^oQUAMvvH}8
z$9I;zyR(!J6{Oq|6|IZ^kyp4NphF`!=GU?g-B%(HZbcOx=JA&@`7Pe#$+~o@>(-Z-
z_*+U(d+KhD^cMIS?Y>vftjsOqkqcYSskZlv1R3Wv@A+lY;l_S%yVstYEtX}MG6Muf
zd^ivCELTvDxV^5v!+KA59|zB?V-~3w6{P<A{J6fBVV!678NnQHm4v7iCBZ(EDfWLj
zT)zKZ720@bPglz>?L|2sn)AQ!J|ioCBuaWp^Nld=xgBc`Olz6!a9o%DkIseER;Ed3
z+K+zAJG)`J##fh3nMNDKu6DflaTgA{xG~d9_Mi9R$5lU){NGI!DXSK0=Gsz!R6IP|
zcj3LWydE+fszJLKm!(&=8eLoZb;7Np4cSHu)~MMi`*gmY#j)VVmK!UkXMCJ=#Vy6?
z)NGf@PIrGz-K}VobMyuOzT2XKOU;h{GYha$ZdmoS!?kF6NP0}%+5a~C?Dp=NF@K@6
z&d#kq!haT@e;_h3-R0XO)<uh3dk$Z0PqcXV?BBx;#zs;X59D6va(_Ib=lWUqH|9%j
z-tcwjUaTx-ZvXZC2kw2#g!1GfjZXwS?Dc&W-fz8xZ_)JWF(+QFJsOu{zK?z4<t5^8
zSwi1L6nU~d*b^l6<Cml3oB6NTWS4Bx?f$K_%WO?1*WKM^M|ZlVn<|}rZI&VxaQ|-F
zwbyaw6aLD{a;>#YU9VJ`dEv(b7Llc!Uj04WI?c`AvEAbh=j+|K|J>Mr+Vk(E-H*3T
zvZ!957NYgeMcSflLFo647J-HK=Ng2h&Xk32WnHLITF~>hTIo&G=Fbef{ca2YS^o5P
ztNZIBm#&z%rvjhJ6<ejR_PpxcSnv3^l<8a3r_x_vcZ>drXJ5t2$@#kZtg58;@xquC
zR=$9q7^w#JKi@fvH9c<sh-zNN_BH#jbd9yI+B3^~0m*%TqdP8#tg7l=*05Lb|7oeo
z?9XN&;?LX>@IU+e_Lj?cSC<z@``0DdKYEgK(Z5c!P)cr3-@&=dt2^qs&iunVY`)!z
z>(B;928L{A1_l8J4#t%H<bb05f@1x|f&x(6P+<|v=(a5e28KD459ZwlZS`aR#@F_~
z_fDF0n?qRO66valjSPj==_yu%THb9(BtFhd&)+02ox8a2^4lxYC+GkB+d1XQ+m+iq
zyA`G0fBI52aVe+OTkpr;Du0F4?Ef`OZQ_Bc340IQyfs+1{ncakot*m)I<O@l64QJ3
zvT;Qur{Gt{><;c}vf5MmO$C-G>eqxy_I{6(&n{8P+w`XUfbc@?jeD}x4wsa~Z2b~<
z`O@oo`gY>6hYxs(%{w5r%x$V?q2H6(%Mptke|&j$@x_D1>otCsUCLNGqkGcI?JN2Z
z23)T`6TI_kLg88a5YZbQL36l2F)oyU#^%%EwQWmQ;ADPJll>KQ%riGNci5)Q+a|zc
zk#uKYu20(yvx!Z+xUEe@%D<h@SK<jOxW0nt;h*Np9J`fIFB~|Snw)QNX+y!uNQq4n
z^<mN+U*z9Bm8njUSm-|ItzpQ5N}lusnpIy^qqy&?{#Oy-kSyu`;OE1(g@(_hH*8%m
zT_DRH{N?bD=|8TB{FCE|7cYvOa+{e&r?UO>p=C`q4`Nwv&tAj-;9JDLZ-2DETz_6!
zWBuyaXQk-08sjum&2#Hru&-;$|5Bf`i;02Z2@eB<I;5phP?VWhl3Em;T2caPYv`l#
z4fT>T^WKCX%)9I$Q0o_b>k`LQ>Em1z9GE&;Ry7?txl!0OFe!SYa%fD~LGSE>+|z|F
z51TALu-@5G^z6{}f)|EDe;L+=&M16*tTa|D)HP&<@LJo=pTC`bbL3vE_J)Vb3%2!$
z#om8l%96Bcm%)ZEt-sSeyB|z4SBrP~7OP>>=BoAhl-@J#Q$HTZaB~0Hw8QAVgOlqm
z`ROnA*-v|+m8|MG|Gew03m^6#vGj?W7n_^;clM3RdDGpC<MXe|{n7mO)miRbxmIv3
zoAGmtcJEc);f-(p?{javx%JkZr$TierzGzQFL>I1c*nb625$mfd=+Irua(|^$7y}y
zs=0i>_DOd#W>-Dovfq|${^8<8)2ksVi#IOR;Fk|y$TBbeIMXgSo7s!P%zN@4$lUs>
zP~`HY^B_ai1MU@vF9=&-C|TQB_TkZiE%#bY4u~&%&7r=~a9fkbg~K=Qmfi2>sM<WU
zWVIXj#jLHnoOx{xIrbz@nj+?Mdt3Xuf*5N>V_mP{cPk?{?+rI(kK1g<bJHmLyp(j&
z{Y#eb6}*gmx3ugNzo6p!*0WS-!l|U7niJZc(elSaHmZ1Bx@%PW>5r@AS;qcfZ>`oV
zt+{?Xt9)CLzy%HF`yNlHaoH<>wp}2@{QcyH6}HzK11s*Fdo1yLHN*UG%WEp{RAufh
zVb-eJ(r1?V_p<o$CEpXI6>V;DP7D04KW*bu>*=dMrm`**ePGP`{_D%S`FF*GLOS;d
z%`F#`aA(M~(bIo<vt_4N&uoS+i<?LPeepS5Htn+F?8i%AbN&8Y*5&_LY)Q$I*Tu)8
zcAMFH-#yLi^ZDuTikd5xTNOMm{_&e-%<o@Y`L4!ybEU{x-Dv4!btksw=v=E0I#_&J
zQzz2yyX&h5cI*Co9-6Tx<!sH_8^?FfQQIu(5H<VF=OU}ag-0UKnYr(0Ow-=Ir8@3N
ze#3izd&~V;3;7nVXJ<B$l>Ch#$=^TwN`S!s)ZJU7cYExu)X~z=*wMX|W8S9OdS@^B
zI!$VQw6ynG=Hg7DO_6Kjj+}pRc6Iu}?)yib`&oZ8?wMAqcDH16X4b+GolB8wGb=w=
zf6kx$J?>V4fPzX~lh<|g3f{(+Epj|-7C7!#os?9O;(J1Da-Dqctivjyg}MA&tE#^g
z`W>>d)7QBk<FUfKeKFTves`hRRGsB7x_5dco7b!_;7iu43J+OZFU@}G?yG-yrdywN
zukWh+e1rYD^}<e@=9%1U9-g}T^tHpkdV7VVuWrVkD-yI{^3C#|umkU+#7m!;B>Yq!
z^$D!|dd)j7=Wxu{DPK#!FzbskPy2F1*zu7{oJYdlExwUkFY?vi_hc#QmDwY7c3R>W
zr|%K@tsAwI;(6FE9CC1wblZ2Mby@c67j5VJLw?7)2s0&Ct}6&OKUOxQrD{rr^}T9A
z=l3O-5|0K5t~Zg|eZ=6X!k6|;vnYx87IsyKzRtNaGxhW=rP*B9Wd5ffYP|l&*Lj=n
z*PE}F-PvAx&bsHaPq1cbep|@g47<!8<{*2WF8)eO)7OVvqNcKWwugGFpT58RujPZn
z{kI>5wmc0Be|h#=|0Ai;2il!yBtIGl88>Dz*{B><yLYMa?~5gIM-0D(bJ|It|07x`
zJ+=3gu7*LT?ArtXmh1PqzCXyP^ka*#&cdJZI?3JN*M&V^&F5ovAf4;p>&<`d_Uost
z;EMdY$N0PjlfmabN~!6>-$K*xNCZ5a?3;7`oXMK!g?cY#^rORX#9i4evuTH_rY~>p
zsmyG1we*$jt(}o_HS)d7darvoO}Sbgb@r!K_r7l1v(E#=qW78i{0(?eqcg2O=rH%w
zDcjb}_|3%oLE`&=lO&&cpRzWl=RA%%x9Fyf!fW4um6o>;dmgDgXLk8M!{)EHw`}zv
zn;+P|{QmYjtV3Ss54t*SCn?bzLlS*__T>PP|4U!TJnuSVF5k|@C7e^p^HDc9Y}O_<
zMbE?^MvHgu%9z$-a`a{Mqx27F!{`6CvpciB@i)W1zAC?KYa*vvN=)v2GiCCe;`8t4
z)%(8Nn|au=VL}C`7P!C{y7bFnLl@iMX&zk<Cz>io`~0filF2tY<mlGeE3e*uTyDwv
z`NyUWLB)y-^<;arOY9}QOQvp|v`GJ%lJA`ldk<Oq?069yxb<M}8}IVf)y47eP2~S+
z{yM8HKj;0FWuJJ`^ft6<hyRab|MdU2oo3kmcZJisf9P#J|8X8emO$^MY2pV8lmyFL
z3f69q`;^4}Y}V2>>yrg4HcD(-a(o_#kord_rq|1IPp#@Ld-2;(X+idr2U9b)@hq0#
z<zA)T6D#;T$)ZJ!jp>xer=xr(vC+=@HFYa)y;qtgz;U{`b@nd9ZHH8@EIxAk+r4gy
zs^ydUwA;EWwqCrU*rwUG(B5nF#X}XhKVB4H_k77Yo#<;KH$yIz{+Z6n{_fW7ww&l!
zm$&*>CcK_gF0p*(vY=P@cvf0nthqCTWl3Fxi2Nt7tZVL^TeWy6@vWZbH*LRvb@744
zj_Z%E<~Y4}-lh3#<&T-JI$$m|)A&*Pk~F5PEENu2zJGk!zb%?vC3yTsyy}k}o6q4#
zazo6twl6raq;w1O@1S^@?t9I-E*~=0Bb+|(k2ogwK6dS~aM>AU4d+GnTu=X9TNi&K
zKy=fm%Jg}KED6slT~D1;d$ao74$cK<)MwwAH#all{NtEQmT}wGZP>da-7=-nJIKoH
zi{<jte4q0{ry)tczh8R&bQaCk@3+kQTq;&8UOnsl;x*f9@+E$6Jn%C@=l_y6@l%?)
zk%r${r5+f*`=5DahRxHfN#}1I+d0Q|v!%ne*?&G4ZEK(0@p(?}ru>GEJDyi>{?dCp
zZ+R_qPKe%n=?8}=Kb6N`;F|<=%Z4&CF!V7oFbE>I>XY+}Qb(Khh@&TZCVA#J8;G>L
zpZZIF!utfC*a@sIu~UMTUA012-_!Z6>UfmxW}fH2+Otu+WU98L-+Mk&efjFeYxeM8
zShgm#%|k$mUu3=O%5RcYTG#5DY?5aRtyv|S@`8nJ?WfP1r<hN<F12k}p2KaSp64E?
z0`8@}NXT_IH2cFaY32dTDZ=Gjr=8ig_}=5TkLKIhWm+GKPD_}XFzq4V3v*+^UXO$6
zDWBJ^mQ%23<5JPcHsIv=_bKYar9YMzw7Hk-yj@v1_t*~aZ4a6)cK0uv)0=ynE%(mx
zf}fmqH*IXbh5ijH-TkQHQ%J!nQ-kQ&f7~=ZK0aC$9e1?n?hf1cUIN>fU;O;LZc|?A
zN}UcHF3(5p8|7og1s6`gZlbbYQ|`Ii<mx@|tyCBu?)*~n!b<b}$Gf!~-Mr;GGv{3B
zKJl^cOnENb0vX}|?K4Ca7pg1?+@#>wu!sNeU-qZ-;iG(vOu7uX&b~qfK7@tkB&+~$
zR0Fo69K{3C1|c^wGcZ7wrGwWYqiaMxo&=%=LSAQM039^}(TjZG1Xgp<52=8d0U{+i
z!R8>FfbNKHw37uu+8}r#H&}BC+BpO02B2^2glGnlsX|}_K(>RI<rgDu@<cZUeQiHT
zKLjt6K{o{wEMOzSi{`Kf3i?uch+!adg)ANe5ew(hU5LK$9Ap#(ua?7OAb2G<Ru_Io
zTgDAB3`DMzM>i1cS+G&)&O%>V2+|M1iK<vl0S6hnDd;nB5EDRTy%E?H#C#mOR`hXl
zkR}KgF^6bHA2mle0d<rIq8CEWuwh_GK^yZy*NxtLhUf&5b@pKU3yA1TqZ^9eCjl7;
z!L?4r7z*y8pc{(b=7SgrB8yy!F%;YaL^l+*G6fk1!+9P^fd{Kv1H4(;KnjEygc$y@
KFfjb`0`UOSKBaj8
literal 0
HcmV?d00001
diff --git a/python/CONTROL.temp b/python/CONTROL.temp
index 2ee9992..c2db90e 100644
--- a/python/CONTROL.temp
+++ b/python/CONTROL.temp
@@ -9,12 +9,12 @@ STREAM OPER
NUMBER OFF
EXPVER 1
GRID 5000
-LEFT -175000
-LOWER -90000
-UPPER 90000
-RIGHT 180000
+LEFT -15000
+LOWER 30000
+UPPER 75000
+RIGHT 45000
LEVEL 60
-LEVELIST 1/to/60
+LEVELIST 55/to/60
RESOL 63
GAUSS 1
ACCURACY 16
@@ -26,8 +26,8 @@ DPDETA 1
SMOOTH 0
FORMAT GRIB1
ADDPAR 186/187/188/235/139/39
-PREFIX EN
-ECSTORAGE 1
+PREFIX EI
+ECSTORAGE 0
ECTRANS 0
ECFSDIR ectmp:/${USER}/econdemand/
MAILFAIL ${USER}
diff --git a/python/Control.py b/python/Control.py
new file mode 100644
index 0000000..8fb6a45
--- /dev/null
+++ b/python/Control.py
@@ -0,0 +1,259 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#************************************************************************
+# TODO AP
+# - write a test class
+# - check documentation
+#************************************************************************
+"""
+@Author: Leopold Haimberger (University of Vienna)
+
+@Date: November 2015
+
+@ChangeHistory:
+ February 2018 - Anne Philipp (University of Vienna):
+ - applied PEP8 style guide
+ - added documentation
+ - applied some minor modifications in programming style/structure
+ - outsource of class Control
+
+@License:
+ (C) Copyright 2015-2018.
+
+ This software is licensed under the terms of the Apache Licence Version 2.0
+ which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+
+@Requirements:
+ A standard python 2.6 or 2.7 installation
+
+@Description:
+ The Control files are the steering part of the FLEXPART extraction
+ software. All necessary parameters needed to retrieve the data fields
+ from the MARS archive for driving FLEXPART are set in a Control file.
+ Some specific parameters like the start and end dates can be overwritten
+ by the command line parameters, but in generel all parameters needed
+ for a complete set of fields for FLEXPART can be set in the Control files.
+
+"""
+# ------------------------------------------------------------------------------
+# MODULES
+# ------------------------------------------------------------------------------
+import os
+import inspect
+import Tools
+# ------------------------------------------------------------------------------
+# CLASS
+# ------------------------------------------------------------------------------
+class Control:
+ '''
+ Class containing the information of the ECMWFDATA control file.
+
+ Contains all the parameters of control files, which are e.g.:
+ DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
+ STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
+ LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
+ OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
+ ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
+ MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR,
+ BASETIME, DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
+
+ For more information about format and content of the parameter
+ see documentation.
+
+ '''
+
+ def __init__(self, filename):
+ '''
+ @Description:
+ Initialises the instance of Control class and defines and
+ assign all controlfile variables. Set default values if
+ parameter was not in CONTROL file.
+
+ @Input:
+ self: instance of Control class
+ Description see class documentation.
+
+ filename: string
+ Name of control file.
+
+ @Return:
+ <nothing>
+ '''
+
+ # read whole CONTROL file
+ with open(filename) as f:
+ fdata = f.read().split('\n')
+
+ # go through every line and store parameter
+ # as class variable
+ for ldata in fdata:
+ data = ldata.split()
+ if len(data) > 1:
+ if 'm_' in data[0].lower():
+ data[0] = data[0][2:]
+ if data[0].lower() == 'class':
+ data[0] = 'marsclass'
+ if data[0].lower() == 'day1':
+ data[0] = 'start_date'
+ if data[0].lower() == 'day2':
+ data[0] = 'end_date'
+ if data[0].lower() == 'addpar':
+ if '/' in data[1]:
+ # remove leading '/' sign from addpar content
+ if data[1][0] == '/':
+ data[1] = data[1][1:]
+ dd = data[1].split('/')
+ data = [data[0]]
+ for d in dd:
+ data.append(d)
+ pass
+ if len(data) == 2:
+ if '$' in data[1]:
+ setattr(self, data[0].lower(), data[1])
+ while '$' in data[1]:
+ i = data[1].index('$')
+ j = data[1].find('{')
+ k = data[1].find('}')
+ var = os.getenv(data[1][j+1:k])
+ if var is not None:
+ data[1] = data[1][:i] + var + data[1][k+1:]
+ else:
+ Tools.myerror(None,
+ 'Could not find variable ' +
+ data[1][j+1:k] +
+ ' while reading ' +
+ filename)
+ setattr(self, data[0].lower() + '_expanded', data[1])
+ else:
+ if data[1].lower() != 'none':
+ setattr(self, data[0].lower(), data[1])
+ else:
+ setattr(self, data[0].lower(), None)
+ elif len(data) > 2:
+ setattr(self, data[0].lower(), (data[1:]))
+ else:
+ pass
+
+ # check a couple of necessary attributes if they contain values
+ # otherwise set default values
+ if not hasattr(self, 'start_date'):
+ self.start_date = None
+ if not hasattr(self, 'end_date'):
+ self.end_date = self.start_date
+ if not hasattr(self, 'accuracy'):
+ self.accuracy = 24
+ if not hasattr(self, 'omega'):
+ self.omega = '0'
+ if not hasattr(self, 'cwc'):
+ self.cwc = '0'
+ if not hasattr(self, 'omegadiff'):
+ self.omegadiff = '0'
+ if not hasattr(self, 'etadiff'):
+ self.etadiff = '0'
+ if not hasattr(self, 'levelist'):
+ if not hasattr(self, 'level'):
+ print('Warning: neither levelist nor level \
+ specified in CONTROL file')
+ else:
+ self.levelist = '1/to/' + self.level
+ else:
+ if 'to' in self.levelist:
+ self.level = self.levelist.split('/')[2]
+ else:
+ self.level = self.levelist.split('/')[-1]
+
+ if not hasattr(self, 'maxstep'):
+ # find out maximum step
+ self.maxstep = 0
+ for s in self.step:
+ if int(s) > self.maxstep:
+ self.maxstep = int(s)
+ else:
+ self.maxstep = int(self.maxstep)
+
+ if not hasattr(self, 'prefix'):
+ self.prefix = 'EN'
+ if not hasattr(self, 'makefile'):
+ self.makefile = None
+ if not hasattr(self, 'basetime'):
+ self.basetime = None
+ if not hasattr(self, 'date_chunk'):
+ self.date_chunk = '3'
+ if not hasattr(self, 'grib2flexpart'):
+ self.grib2flexpart = '0'
+
+ # script directory
+ self.ecmwfdatadir = os.path.dirname(os.path.abspath(inspect.getfile(
+ inspect.currentframe()))) + '/../'
+ # Fortran source directory
+ self.exedir = self.ecmwfdatadir + 'src/'
+
+ # FLEXPART directory
+ if not hasattr(self, 'flexpart_root_scripts'):
+ self.flexpart_root_scripts = self.ecmwfdatadir
+
+ return
+
+ def __str__(self):
+ '''
+ @Description:
+ Prepares a single string with all the comma seperated Control
+ class attributes including their values.
+
+ Example:
+ {'kids': 0, 'name': 'Dog', 'color': 'Spotted',
+ 'age': 10, 'legs': 2, 'smell': 'Alot'}
+
+ @Input:
+ self: instance of Control class
+ Description see class documentation.
+
+ @Return:
+ string of Control class attributes with their values
+ '''
+
+ attrs = vars(self)
+
+ return ', '.join("%s: %s" % item for item in attrs.items())
+
+ def tolist(self):
+ '''
+ @Description:
+ Just generates a list of strings containing the attributes and
+ assigned values except the attributes "_expanded", "exedir",
+ "ecmwfdatadir" and "flexpart_root_scripts".
+
+ @Input:
+ self: instance of Control class
+ Description see class documentation.
+
+ @Return:
+ l: list
+ A sorted list of the all Control class attributes with
+ their values except the attributes "_expanded", "exedir",
+ "ecmwfdatadir" and "flexpart_root_scripts".
+ '''
+
+ attrs = vars(self)
+ l = list()
+
+ for item in attrs.items():
+ if '_expanded' in item[0]:
+ pass
+ elif 'exedir' in item[0]:
+ pass
+ elif 'flexpart_root_scripts' in item[0]:
+ pass
+ elif 'ecmwfdatadir' in item[0]:
+ pass
+ else:
+ if type(item[1]) is list:
+ stot = ''
+ for s in item[1]:
+ stot += s + ' '
+
+ l.append("%s %s" % (item[0], stot))
+ else:
+ l.append("%s %s" % item)
+
+ return sorted(l)
diff --git a/python/Disagg.py b/python/Disagg.py
new file mode 100644
index 0000000..cf2f4c9
--- /dev/null
+++ b/python/Disagg.py
@@ -0,0 +1,156 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#************************************************************************
+# TODO AP
+# - make a class out of this ???
+# - write a test class
+#************************************************************************
+"""
+@Author: Anne Philipp (University of Vienna)
+
+@Date: March 2018
+
+@ChangeHistory:
+ November 2015 - Leopold Haimberger (University of Vienna):
+ - integrated methods dapoly and darain from Fortran to Python
+
+ April 2018 - Anne Philipp (University of Vienna):
+ - applied PEP8 style guide
+ - added structured documentation
+ - outsourced the disaggregation functions dapoly and darain
+ to a new module named Disagg
+
+@License:
+ (C) Copyright 2015-2018.
+
+ This software is licensed under the terms of the Apache Licence Version 2.0
+ which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+
+@Requirements:
+ A standard python 2.6 or 2.7 installation
+
+@Description:
+ Further documentation may be obtained from www.flexpart.eu.
+
+ Functionality provided:
+ Disaggregation of deaccumulated flux data from an ECMWF model FG field.
+ Initially the flux data to be concerned are:
+ - large-scale precipitation
+ - convective precipitation
+ - surface sensible heat flux
+ - surface solar radiation
+ - u stress
+ - v stress
+
+ Different versions of disaggregation is provided for rainfall
+ data (darain, modified linear) and the surface fluxes and
+ stress data (dapoly, cubic polynomial).
+
+"""
+# ------------------------------------------------------------------------------
+# MODULES
+# ------------------------------------------------------------------------------
+
+# ------------------------------------------------------------------------------
+# FUNCTIONS
+# ------------------------------------------------------------------------------
+def dapoly(alist):
+ '''
+ @Author: P. JAMES
+
+ @Date: 2000-03-29
+
+ @ChangeHistory:
+ June 2003 - A. BECK (2003-06-01)
+ adaptaions
+ November 2015 - Leopold Haimberger (University of Vienna)
+ migration from Fortran to Python
+
+ @Description:
+ Interpolation of deaccumulated fluxes of an ECMWF model FG field
+ using a cubic polynomial solution which conserves the integrals
+ of the fluxes within each timespan.
+ Disaggregation is done for 4 accumluated timespans which generates
+ a new, disaggregated value which is output at the central point
+ of the 4 accumulation timespans. This new point is used for linear
+ interpolation of the complete timeseries afterwards.
+
+ @Input:
+ alist: list of size 4, array(2D), type=float
+ List of 4 timespans as 2-dimensional, horizontal fields.
+ E.g. [[array_t1], [array_t2], [array_t3], [array_t4]]
+
+ @Return:
+ nfield: array(2D), type=float
+ New field which replaces the field at the second position
+ of the accumulation timespans.
+
+ '''
+ pya = (alist[3] - alist[0] + 3. * (alist[1] - alist[2])) / 6.
+ pyb = (alist[2] + alist[0]) / 2. - alist[1] - 9. * pya / 2.
+ pyc = alist[1] - alist[0] - 7. * pya / 2. - 2. * pyb
+ pyd = alist[0] - pya / 4. - pyb / 3. - pyc / 2.
+ nfield = 8. * pya + 4. * pyb + 2. * pyc + pyd
+
+ return nfield
+
+
+def darain(alist):
+ '''
+ @Author: P. JAMES
+
+ @Date: 2000-03-29
+
+ @ChangeHistory:
+ June 2003 - A. BECK (2003-06-01)
+ adaptaions
+ November 2015 - Leopold Haimberger (University of Vienna)
+ migration from Fortran to Python
+
+ @Description:
+ Interpolation of deaccumulated fluxes of an ECMWF model FG rainfall
+ field using a modified linear solution which conserves the integrals
+ of the fluxes within each timespan.
+ Disaggregation is done for 4 accumluated timespans which generates
+ a new, disaggregated value which is output at the central point
+ of the 4 accumulation timespans. This new point is used for linear
+ interpolation of the complete timeseries afterwards.
+
+ @Input:
+ alist: list of size 4, array(2D), type=float
+ List of 4 timespans as 2-dimensional, horizontal fields.
+ E.g. [[array_t1], [array_t2], [array_t3], [array_t4]]
+
+ @Return:
+ nfield: array(2D), type=float
+ New field which replaces the field at the second position
+ of the accumulation timespans.
+ '''
+ xa = alist[0]
+ xb = alist[1]
+ xc = alist[2]
+ xd = alist[3]
+ xa[xa < 0.] = 0.
+ xb[xb < 0.] = 0.
+ xc[xc < 0.] = 0.
+ xd[xd < 0.] = 0.
+
+ xac = 0.5 * xb
+ mask = xa + xc > 0.
+ xac[mask] = xb[mask] * xc[mask] / (xa[mask] + xc[mask])
+ xbd = 0.5 * xc
+ mask = xb + xd > 0.
+ xbd[mask] = xb[mask] * xc[mask] / (xb[mask] + xd[mask])
+ nfield = xac + xbd
+
+ return nfield
+
+
+
+
+
+
+
+
+
+
diff --git a/python/ECFlexpart.py b/python/ECFlexpart.py
new file mode 100644
index 0000000..8e1e57f
--- /dev/null
+++ b/python/ECFlexpart.py
@@ -0,0 +1,1253 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#************************************************************************
+# TODO AP
+#AP
+# - specifiy file header documentation
+# - apply classtests
+# - add references to ECMWF specific software packages
+#************************************************************************
+"""
+@Author: Anne Fouilloux (University of Oslo)
+
+@Date: October 2014
+
+@ChangeHistory:
+ November 2015 - Leopold Haimberger (University of Vienna):
+ - extended with Class Control
+ - removed functions mkdir_p, daterange, years_between, months_between
+ - added functions darain, dapoly, toparamId, init128, normalexit,
+ myerror, cleanup, install_args_and_control,
+ interpret_args_and_control,
+ - removed function __del__ in class EIFLexpart
+ - added the following functions in EIFlexpart:
+ - create_namelist
+ - process_output
+ - deacc_fluxes
+ - modified existing EIFlexpart - functions for the use in
+ flex_extract
+ - retrieve also longer term forecasts, not only analyses and
+ short term forecast data
+ - added conversion into GRIB2
+ - added conversion into .fp format for faster execution of FLEXPART
+ (see https://www.flexpart.eu/wiki/FpCtbtoWo4FpFormat)
+
+ February 2018 - Anne Philipp (University of Vienna):
+ - applied PEP8 style guide
+ - added documentation
+ - outsourced class Control
+ - outsourced class MarsRetrieval
+ - changed class name from EIFlexpart to ECFlexpart
+ - applied minor code changes (style)
+
+@License:
+ (C) Copyright 2014-2018.
+
+ This software is licensed under the terms of the Apache Licence Version 2.0
+ which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+
+@Requirements:
+ - A standard python 2.6 or 2.7 installation
+ - dateutils
+ - matplotlib (optional, for debugging)
+ - ECMWF specific packages, all available from https://software.ecmwf.int/
+ ECMWF WebMARS, gribAPI with python enabled, emoslib and
+ ecaccess web toolkit
+
+@Description:
+ Further documentation may be obtained from www.flexpart.eu.
+
+ Functionality provided:
+ Prepare input 3D-wind fields in hybrid coordinates +
+ surface fields for FLEXPART runs
+"""
+# ------------------------------------------------------------------------------
+# MODULES
+# ------------------------------------------------------------------------------
+import subprocess
+from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter
+import traceback
+import shutil
+import os
+import errno
+import sys
+import inspect
+import glob
+import datetime
+from string import atoi
+from numpy import *
+ecapi = True
+try:
+ import ecmwfapi
+except ImportError:
+ ecapi = False
+from gribapi import *
+from GribTools import GribTools
+from Tools import init128, toparamId, silentremove, product
+from Control import Control
+from MARSretrieval import MARSretrieval
+import Disagg
+# ------------------------------------------------------------------------------
+# CLASS
+# ------------------------------------------------------------------------------
+class ECFlexpart:
+ '''
+ Class to retrieve ECMWF data specific for running FLEXPART.
+ '''
+ # --------------------------------------------------------------------------
+ # CLASS FUNCTIONS
+ # --------------------------------------------------------------------------
+ def __init__(self, c, fluxes=False): #done/ verstehen
+ '''
+ @Description:
+ Creates an object/instance of ECFlexpart with the
+ associated settings of its attributes for the retrieval.
+
+ @Input:
+ self: instance of ECFlexpart
+ The current object of the class.
+
+ c: instance of class Control
+ Contains all the parameters of control files, which are e.g.:
+ DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
+ STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
+ LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
+ OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
+ ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
+ MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
+ DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
+
+ For more information about format and content of the parameter
+ see documentation.
+
+ fluxes: boolean, optional
+ Decides if a the flux parameter settings are stored or
+ the rest of the parameter list.
+ Default value is False.
+
+ @Return:
+ <nothing>
+ '''
+
+ # different mars types for retrieving reanalysis data for flexpart
+ self.types = dict()
+ try:
+ if c.maxstep > len(c.type): # Pure forecast mode
+ c.type = [c.type[1]]
+ c.step = ['{:0>3}'.format(int(c.step[0]))]
+ c.time = [c.time[0]]
+ for i in range(1, c.maxstep + 1):
+ c.type.append(c.type[0])
+ c.step.append('{:0>3}'.format(i))
+ c.time.append(c.time[0])
+ except:
+ pass
+
+ self.inputdir = c.inputdir
+ self.basetime = c.basetime
+ self.dtime = c.dtime
+ #self.mars = {}
+ i = 0
+ #j = 0
+ if fluxes is True and c.maxstep < 24:
+ # no forecast beyond one day is needed!
+ # Thus, prepare flux data manually as usual
+ # with only FC fields with start times at 00/12
+ # (but without 00/12 fields since these are
+ # the initialisation times of the flux fields
+ # and therefore are zero all the time)
+ self.types['FC'] = {'times': '00/12',
+ 'steps': '{}/to/12/by/{}'.format(c.dtime,
+ c.dtime)}
+ #i = 1
+ #for k in [0, 12]:
+ # for j in range(int(c.dtime), 13, int(c.dtime)):
+ # self.mars['{:0>3}'.format(i * int(c.dtime))] = \
+ # [c.type[1], '{:0>2}'.format(k), '{:0>3}'.format(j)]
+ # i += 1
+ else:
+ for ty, st, ti in zip(c.type, c.step, c.time):
+ btlist = range(24)
+ if c.basetime == '12':
+ btlist = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
+ if c.basetime == '00':
+ btlist = [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 0]
+
+ if mod(i, int(c.dtime)) == 0 and \
+ (c.maxstep > 24 or i in btlist):
+
+ if ty not in self.types.keys():
+ self.types[ty] = {'times': '', 'steps': ''}
+
+ if ti not in self.types[ty]['times']:
+ if len(self.types[ty]['times']) > 0:
+ self.types[ty]['times'] += '/'
+ self.types[ty]['times'] += ti
+
+ if st not in self.types[ty]['steps']:
+ if len(self.types[ty]['steps']) > 0:
+ self.types[ty]['steps'] += '/'
+ self.types[ty]['steps'] += st
+
+ #self.mars['{:0>3}'.format(j)] = [ty,
+ # '{:0>2}'.format(int(ti)),
+ # '{:0>3}'.format(int(st))]
+ #j += int(c.dtime)
+
+ i += 1
+ print 'EC init: ', self.types #AP
+ # Different grids need different retrievals
+ # SH = Spherical Harmonics, GG = Gaussian Grid,
+ # OG = Output Grid, ML = MultiLevel, SL = SingleLevel
+ self.params = {'SH__ML': '', 'SH__SL': '',
+ 'GG__ML': '', 'GG__SL': '',
+ 'OG__ML': '', 'OG__SL': '',
+ 'OG_OROLSM_SL': '', 'OG_acc_SL': ''}
+ self.marsclass = c.marsclass
+ self.stream = c.stream
+ self.number = c.number
+ self.resol = c.resol
+ self.accuracy = c.accuracy
+ self.level = c.level
+ try:
+ self.levelist = c.levelist
+ except:
+ self.levelist = '1/to/' + c.level
+
+ # for gaussian grid retrieval
+ self.glevelist = '1/to/' + c.level
+
+ try:
+ self.gaussian = c.gaussian
+ except:
+ self.gaussian = ''
+
+ try:
+ self.expver = c.expver
+ except:
+ self.expver = '1'
+
+ try:
+ self.number = c.number
+ except:
+ self.number = '0'
+
+ if 'N' in c.grid: # Gaussian output grid
+ self.grid = c.grid
+ self.area = 'G'
+ else:
+ self.grid = '{}/{}'.format(int(c.grid) / 1000., int(c.grid) / 1000.)
+ self.area = '{}/{}/{}/{}'.format(int(c.upper) / 1000.,
+ int(c.left) / 1000.,
+ int(c.lower) / 1000.,
+ int(c.right) / 1000.)
+
+ self.outputfilelist = []
+
+
+ # Now comes the nasty part that deals with the different
+ # scenarios we have:
+ # 1) Calculation of etadot on
+ # a) Gaussian grid
+ # b) Output grid
+ # c) Output grid using parameter 77 retrieved from MARS
+ # 3) Calculation/Retrieval of omega
+ # 4) Download also data for WRF
+
+ if fluxes is False:
+ self.params['SH__SL'] = ['LNSP', 'ML', '1', 'OFF']
+ # "SD/MSL/TCC/10U/10V/2T/2D/129/172"
+ self.params['OG__SL'] = ["141/151/164/165/166/167/168/129/172", \
+ 'SFC', '1', self.grid]
+ if len(c.addpar) > 0:
+ if c.addpar[0] == '/':
+ c.addpar = c.addpar[1:]
+ self.params['OG__SL'][0] += '/' + '/'.join(c.addpar)
+
+ self.params['OG_OROLSM__SL'] = ["160/27/28/173", \
+ 'SFC', '1', self.grid]
+
+ self.params['OG__ML'] = ['T/Q', 'ML', self.levelist, self.grid]
+
+ if c.gauss == '0' and c.eta == '1':
+ # the simplest case
+ self.params['OG__ML'][0] += '/U/V/77'
+ elif c.gauss == '0' and c.eta == '0':
+#AP then remove?!?!?!? # this is not recommended (inaccurate)
+ self.params['OG__ML'][0] += '/U/V'
+ elif c.gauss == '1' and c.eta == '0':
+ # this is needed for data before 2008, or for reanalysis data
+ self.params['GG__SL'] = ['Q', 'ML', '1', \
+ '{}'.format((int(self.resol) + 1) / 2)]
+ self.params['SH__ML'] = ['U/V/D', 'ML', self.glevelist, 'OFF']
+ else:
+ print('Warning: This is a very costly parameter combination, \
+ use only for debugging!')
+ self.params['GG__SL'] = ['Q', 'ML', '1', \
+ '{}'.format((int(self.resol) + 1) / 2)]
+ self.params['GG__ML'] = ['U/V/D/77', 'ML', self.glevelist, \
+ '{}'.format((int(self.resol) + 1) / 2)]
+
+ if c.omega == '1':
+ self.params['OG__ML'][0] += '/W'
+
+ try:
+ # add cloud water content if necessary
+ if c.cwc == '1':
+ self.params['OG__ML'][0] += '/CLWC/CIWC'
+ except:
+ pass
+
+ try:
+ # add vorticity and geopotential height for WRF if necessary
+ if c.wrf == '1':
+ self.params['OG__ML'][0] += '/Z/VO'
+ if '/D' not in self.params['OG__ML'][0]:
+ self.params['OG__ML'][0] += '/D'
+ #wrf_sfc = 'sp/msl/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/stl1/ /
+ # stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4'.upper()
+ wrf_sfc = '134/235/167/165/166/168/129/172/34/31/141/ \
+ 139/170/183/236/39/40/41/42'.upper()
+ lwrt_sfc = wrf_sfc.split('/')
+ for par in lwrt_sfc:
+ if par not in self.params['OG__SL'][0]:
+ self.params['OG__SL'][0] += '/' + par
+ except:
+ pass
+ else:
+ self.params['OG_acc_SL'] = ["LSP/CP/SSHF/EWSS/NSSS/SSR", \
+ 'SFC', '1', self.grid]
+
+ # if needed, add additional WRF specific parameters here
+
+ return
+
+
+ def write_namelist(self, c, filename): #done
+ '''
+ @Description:
+ Creates a namelist file in the temporary directory and writes
+ the following values to it: maxl, maxb, mlevel,
+ mlevelist, mnauf, metapar, rlo0, rlo1, rla0, rla1,
+ momega, momegadiff, mgauss, msmooth, meta, metadiff, mdpdeta
+
+ @Input:
+ self: instance of ECFlexpart
+ The current object of the class.
+
+ c: instance of class Control
+ Contains all the parameters of control files, which are e.g.:
+ DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
+ STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
+ LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
+ OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
+ ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
+ MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
+ DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
+
+ For more information about format and content of the parameter
+ see documentation.
+
+ filename: string
+ Name of the namelist file.
+
+ @Return:
+ <nothing>
+ '''
+
+ self.inputdir = c.inputdir
+ area = asarray(self.area.split('/')).astype(float)
+ grid = asarray(self.grid.split('/')).astype(float)
+
+ if area[1] > area[3]:
+ area[1] -= 360
+ zyk = abs((area[3] - area[1] - 360.) + grid[1]) < 1.e-6
+ maxl = int((area[3] - area[1]) / grid[1]) + 1
+ maxb = int((area[0] - area[2]) / grid[0]) + 1
+
+ with open(self.inputdir + '/' + filename, 'w') as f:
+ f.write('&NAMGEN\n')
+ f.write(',\n '.join(['maxl = ' + str(maxl), 'maxb = ' + str(maxb),
+ 'mlevel = ' + self.level,
+ 'mlevelist = ' + '"' + self.levelist + '"',
+ 'mnauf = ' + self.resol, 'metapar = ' + '77',
+ 'rlo0 = ' + str(area[1]), 'rlo1 = ' + str(area[3]),
+ 'rla0 = ' + str(area[2]), 'rla1 = ' + str(area[0]),
+ 'momega = ' + c.omega, 'momegadiff = ' + c.omegadiff,
+ 'mgauss = ' + c.gauss, 'msmooth = ' + c.smooth,
+ 'meta = ' + c.eta, 'metadiff = ' + c.etadiff,
+ 'mdpdeta = ' + c.dpdeta]))
+
+ f.write('\n/\n')
+
+ return
+
+ def retrieve(self, server, dates, times, inputdir=''):
+ '''
+ @Description:
+
+
+ @Input:
+ self: instance of ECFlexpart
+
+ server: instance of ECMWFService
+
+ dates:
+
+ times:
+
+ inputdir: string, optional
+ Default string is empty ('').
+
+ @Return:
+ <nothing>
+ '''
+ self.dates = dates
+ self.server = server
+
+ if inputdir == "":
+ self.inputdir = '.'
+ else:
+ self.inputdir = inputdir
+
+ # Retrieve Q not for using Q but as a template for a reduced gaussian
+ # grid one date and time is enough
+ # Take analysis at 00
+ qdate = self.dates
+ idx = qdate.find("/")
+ if (idx > 0):
+ qdate = self.dates[:idx]
+
+ #QG = MARSretrieval(self.server, marsclass = self.marsclass, stream = self.stream, type = "an", levtype = "ML", levelist = "1",
+ #gaussian = "reduced",grid = '{}'.format((int(self.resol)+1)/2), resol = self.resol,accuracy = self.accuracy,target = self.inputdir+"/"+"QG.grb",
+ #date = qdate, time = "00",expver = self.expver, param = "133.128")
+ #QG.displayInfo()
+ #QG.dataRetrieve()
+
+ oro = False
+ for ftype in self.types:
+ for pk, pv in self.params.iteritems():
+ if isinstance(pv, str):
+ continue
+ mftype = '' + ftype
+ mftime = self.types[ftype]['times']
+ mfstep = self.types[ftype]['steps']
+ mfdate = self.dates
+ mfstream = self.stream
+ mftarget = self.inputdir + "/" + ftype + pk + '.' + \
+ self.dates.split('/')[0] + '.' + str(os.getppid()) +\
+ '.' + str(os.getpid()) + ".grb"
+ if pk == 'OG__SL':
+ pass
+ if pk == 'OG_OROLSM__SL':
+ if oro is False:
+ mfstream = 'OPER'
+ mftype = 'AN'
+ mftime = '00'
+ mfstep = '000'
+ mfdate = self.dates.split('/')[0]
+ mftarget = self.inputdir + "/" + pk + '.' + mfdate + \
+ '.' + str(os.getppid()) + '.' + \
+ str(os.getpid()) + ".grb"
+ oro = True
+ else:
+ continue
+
+ if pk == 'GG__SL' and pv[0] == 'Q':
+ area = ""
+ gaussian = 'reduced'
+ else:
+ area = self.area
+ gaussian = self.gaussian
+
+ if self.basetime is None:
+ MR = MARSretrieval(self.server,
+ marsclass=self.marsclass, stream=mfstream,
+ type=mftype, levtype=pv[1], levelist=pv[2],
+ resol=self.resol, gaussian=gaussian,
+ accuracy=self.accuracy, grid=pv[3],
+ target=mftarget, area=area, date=mfdate,
+ time=mftime, number=self.number, step=mfstep,
+ expver=self.expver, param=pv[0])
+
+ MR.displayInfo()
+ MR.dataRetrieve()
+ # The whole else section is only necessary for operational scripts.
+ # It could be removed
+ else:
+ # check if mars job requests fields beyond basetime.
+ # If yes eliminate those fields since they may not
+ # be accessible with user's credentials
+ sm1 = -1
+ if 'by' in mfstep:
+ sm1 = 2
+ tm1 = -1
+ if 'by' in mftime:
+ tm1 = 2
+ maxtime = datetime.datetime.strptime(
+ mfdate.split('/')[-1] + mftime.split('/')[tm1],
+ '%Y%m%d%H') + datetime.timedelta(
+ hours=int(mfstep.split('/')[sm1]))
+
+ elimit = datetime.datetime.strptime(
+ mfdate.split('/')[-1] +
+ self.basetime, '%Y%m%d%H')
+
+ if self.basetime == '12':
+ if 'acc' in pk:
+
+ # Strategy: if maxtime-elimit> = 24h reduce date by 1,
+ # if 12h< = maxtime-elimit<12h reduce time for last date
+ # if maxtime-elimit<12h reduce step for last time
+ # A split of the MARS job into 2 is likely necessary.
+ maxtime = elimit-datetime.timedelta(hours=24)
+ mfdate = '/'.join(('/'.join(mfdate.split('/')[:-1]),
+ datetime.datetime.strftime(
+ maxtime, '%Y%m%d')))
+
+ MR = MARSretrieval(self.server,
+ marsclass=self.marsclass,
+ stream=self.stream, type=mftype,
+ levtype=pv[1], levelist=pv[2],
+ resol=self.resol, gaussian=gaussian,
+ accuracy=self.accuracy, grid=pv[3],
+ target=mftarget, area=area,
+ date=mfdate, time=mftime,
+ number=self.number, step=mfstep,
+ expver=self.expver, param=pv[0])
+
+ MR.displayInfo()
+ MR.dataRetrieve()
+
+ maxtime = elimit - datetime.timedelta(hours=12)
+ mfdate = datetime.datetime.strftime(maxtime,
+ '%Y%m%d')
+ mftime = '00'
+ mftarget = self.inputdir + "/" + ftype + pk + \
+ '.' + mfdate + '.' + str(os.getppid()) +\
+ '.' + str(os.getpid()) + ".grb"
+
+ MR = MARSretrieval(self.server,
+ marsclass=self.marsclass,
+ stream=self.stream, type=mftype,
+ levtype=pv[1], levelist=pv[2],
+ resol=self.resol, gaussian=gaussian,
+ accuracy=self.accuracy, grid=pv[3],
+ target=mftarget, area=area,
+ date=mfdate, time=mftime,
+ number=self.number, step=mfstep,
+ expver=self.expver, param=pv[0])
+
+ MR.displayInfo()
+ MR.dataRetrieve()
+ else:
+ MR = MARSretrieval(self.server,
+ marsclass=self.marsclass,
+ stream=self.stream, type=mftype,
+ levtype=pv[1], levelist=pv[2],
+ resol=self.resol, gaussian=gaussian,
+ accuracy=self.accuracy, grid=pv[3],
+ target=mftarget, area=area,
+ date=mfdate, time=mftime,
+ number=self.number, step=mfstep,
+ expver=self.expver, param=pv[0])
+
+ MR.displayInfo()
+ MR.dataRetrieve()
+ else:
+ maxtime = elimit - datetime.timedelta(hours=24)
+ mfdate = datetime.datetime.strftime(maxtime,'%Y%m%d')
+
+ mftimesave = ''.join(mftime)
+
+ if '/' in mftime:
+ times = mftime.split('/')
+ while ((int(times[0]) +
+ int(mfstep.split('/')[0]) <= 12) and
+ (pk != 'OG_OROLSM__SL') and 'acc' not in pk):
+ times = times[1:]
+ if len(times) > 1:
+ mftime = '/'.join(times)
+ else:
+ mftime = times[0]
+
+ MR = MARSretrieval(self.server,
+ marsclass=self.marsclass,
+ stream=self.stream, type=mftype,
+ levtype=pv[1], levelist=pv[2],
+ resol=self.resol, gaussian=gaussian,
+ accuracy=self.accuracy, grid=pv[3],
+ target=mftarget, area=area,
+ date=mfdate, time=mftime,
+ number=self.number, step=mfstep,
+ expver=self.expver, param=pv[0])
+
+ MR.displayInfo()
+ MR.dataRetrieve()
+
+ if (int(mftimesave.split('/')[0]) == 0 and
+ int(mfstep.split('/')[0]) == 0 and
+ pk != 'OG_OROLSM__SL'):
+ mfdate = datetime.datetime.strftime(elimit,'%Y%m%d')
+ mftime = '00'
+ mfstep = '000'
+ mftarget = self.inputdir + "/" + ftype + pk + \
+ '.' + mfdate + '.' + str(os.getppid()) +\
+ '.' + str(os.getpid()) + ".grb"
+
+ MR = MARSretrieval(self.server,
+ marsclass=self.marsclass,
+ stream=self.stream, type=mftype,
+ levtype=pv[1], levelist=pv[2],
+ resol=self.resol, gaussian=gaussian,
+ accuracy=self.accuracy, grid=pv[3],
+ target=mftarget, area=area,
+ date=mfdate, time=mftime,
+ number=self.number, step=mfstep,
+ expver=self.expver, param=pv[0])
+
+ MR.displayInfo()
+ MR.dataRetrieve()
+
+ print("MARS retrieve done... ")
+
+ return
+
+
+ def process_output(self, c): #done
+ '''
+ @Description:
+ The grib files are postprocessed depending on selection in
+ control file. The resulting files are moved to the output
+ directory if its not equla to the input directory.
+ The following modifications might be done if
+ properly switched in control file:
+ GRIB2 - Conversion to GRIB2
+ ECTRANS - Transfer of files to gateway server
+ ECSTORAGE - Storage at ECMWF server
+ GRIB2FLEXPART - Conversion of GRIB files to FLEXPART binary format
+
+ @Input:
+ self: instance of ECFlexpart
+ The current object of the class.
+
+ c: instance of class Control
+ Contains all the parameters of control files, which are e.g.:
+ DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
+ STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
+ LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
+ OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
+ ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
+ MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
+ DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
+
+ For more information about format and content of the parameter
+ see documentation.
+
+ @Return:
+ <nothing>
+
+ '''
+
+ print('Postprocessing:\n Format: {}\n'.format(c.format))
+
+ if c.ecapi is False:
+ print('ecstorage: {}\n ecfsdir: {}\n'.
+ format(c.ecstorage, c.ecfsdir))
+ if not hasattr(c, 'gateway'):
+ c.gateway = os.getenv('GATEWAY')
+ if not hasattr(c, 'destination'):
+ c.destination = os.getenv('DESTINATION')
+ print('ectrans: {}\n gateway: {}\n destination: {}\n '
+ .format(c.ectrans, c.gateway, c.destination))
+
+ print('Output filelist: \n')
+ print(self.outputfilelist)
+
+ if c.format.lower() == 'grib2':
+ for ofile in self.outputfilelist:
+ p = subprocess.check_call(['grib_set', '-s', 'edition=2, \
+ productDefinitionTemplateNumber=8',
+ ofile, ofile + '_2'])
+ p = subprocess.check_call(['mv', ofile + '_2', ofile])
+
+ if int(c.ectrans) == 1 and c.ecapi is False:
+ for ofile in self.outputfilelist:
+ p = subprocess.check_call(['ectrans', '-overwrite', '-gateway',
+ c.gateway, '-remote', c.destination,
+ '-source', ofile])
+ print('ectrans:', p)
+
+ if int(c.ecstorage) == 1 and c.ecapi is False:
+ for ofile in self.outputfilelist:
+ p = subprocess.check_call(['ecp', '-o', ofile,
+ os.path.expandvars(c.ecfsdir)])
+
+ if c.outputdir != c.inputdir:
+ for ofile in self.outputfilelist:
+ p = subprocess.check_call(['mv', ofile, c.outputdir])
+
+ # prepare environment for the grib2flexpart run
+ # to convert grib to flexpart binary
+ if c.grib2flexpart == '1':
+
+ # generate AVAILABLE file
+ # Example of AVAILABLE file data
+ # 20131107 000000 EN13110700 ON DISC
+ clist = []
+ for ofile in self.outputfilelist:
+ fname = ofile.split('/')
+ if '.' in fname[-1]:
+ l = fname[-1].split('.')
+ timestamp = datetime.datetime.strptime(l[0][-6:] + l[1],
+ '%y%m%d%H')
+ timestamp += datetime.timedelta(hours=int(l[2]))
+ cdate = datetime.datetime.strftime(timestamp, '%Y%m%d')
+ chms = datetime.datetime.strftime(timestamp, '%H%M%S')
+ else:
+ cdate = '20' + fname[-1][-8:-2]
+ chms = fname[-1][-2:] + '0000'
+ clist.append(cdate + ' ' + chms + ' '*6 +
+ fname[-1] + ' '*14 + 'ON DISC')
+ clist.sort()
+ with open(c.outputdir + '/' + 'AVAILABLE', 'w') as f:
+ f.write('\n'.join(clist) + '\n')
+
+ # generate pathnames file
+ pwd = os.path.abspath(c.outputdir)
+ with open(pwd + '/pathnames','w') as f:
+ f.write(pwd + '/Options/\n')
+ f.write(pwd + '/\n')
+ f.write(pwd + '/\n')
+ f.write(pwd + '/AVAILABLE\n')
+ f.write(' = == = == = == = == = == == = \n')
+
+ # create Options dir if necessary
+ if not os.path.exists(pwd + '/Options'):
+ os.makedirs(pwd+'/Options')
+
+ # read template COMMAND file
+ with open(os.path.expandvars(
+ os.path.expanduser(c.flexpart_root_scripts)) +
+ '/../Options/COMMAND', 'r') as f:
+ lflist = f.read().split('\n')
+
+ # find index of list where to put in the
+ # date and time information
+ # usually after the LDIRECT parameter
+ i = 0
+ for l in lflist:
+ if 'LDIRECT' in l.upper():
+ break
+ i += 1
+
+# clist.sort()
+ # insert the date and time information of run star and end
+ # into the list of lines of COMMAND file
+ lflist = lflist[:i+1] + \
+ [clist[0][:16], clist[-1][:16]] + \
+ lflist[i+3:]
+
+ # write the new COMMAND file
+ with open(pwd + '/Options/COMMAND', 'w') as g:
+ g.write('\n'.join(lflist) + '\n')
+
+ # change to outputdir and start the
+ # grib2flexpart run
+ # afterwards switch back to the working dir
+ os.chdir(c.outputdir)
+ p = subprocess.check_call([os.path.expandvars(
+ os.path.expanduser(c.flexpart_root_scripts)) +
+ '/../FLEXPART_PROGRAM/grib2flexpart',
+ 'useAvailable', '.'])
+ os.chdir(pwd)
+
+ return
+
+ def create(self, inputfiles, c): #done
+ '''
+ @Description:
+ This method is based on the ECMWF example index.py
+ https://software.ecmwf.int/wiki/display/GRIB/index.py
+
+ An index file will be created which depends on the combination
+ of "date", "time" and "stepRange" values. This is used to iterate
+ over all messages in the grib files passed through the parameter
+ "inputfiles" to seperate specific parameters into fort.* files.
+ Afterwards the FORTRAN program Convert2 is called to convert
+ the data fields all to the same grid and put them in one file
+ per day.
+
+ @Input:
+ self: instance of ECFlexpart
+ The current object of the class.
+
+ inputfiles: instance of UIOFiles
+ Contains a list of files.
+
+ c: instance of class Control
+ Contains all the parameters of control files, which are e.g.:
+ DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
+ STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
+ LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
+ OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
+ ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
+ MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
+ DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
+
+ For more information about format and content of the parameter
+ see documentation.
+
+ @Return:
+ <nothing>
+ '''
+
+ table128 = init128(c.ecmwfdatadir +
+ '/grib_templates/ecmwf_grib1_table_128')
+ wrfpars = toparamId('sp/mslp/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/\
+ stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4',
+ table128)
+
+ index_keys = ["date", "time", "step"]
+ indexfile = c.inputdir + "/date_time_stepRange.idx"
+ silentremove(indexfile)
+ grib = GribTools(inputfiles.files)
+ # creates new index file
+ iid = grib.index(index_keys=index_keys, index_file=indexfile)
+
+ # read values of index keys
+ index_vals = []
+ for key in index_keys:
+ index_vals.append(grib_index_get(iid, key))
+ print(index_vals[-1])
+ # index_vals looks for example like:
+ # index_vals[0]: ('20171106', '20171107', '20171108') ; date
+ # index_vals[1]: ('0', '1200', '1800', '600') ; time
+ # index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange
+
+ # delete old fort.* files and open them newly
+ fdict = {'10':None, '11':None, '12':None, '13':None, '16':None,
+ '17':None, '19':None, '21':None, '22':None, '20':None}
+ #for f in fdict.keys():
+ # silentremove(c.inputdir + "/fort." + f)
+
+
+ for prod in product(*index_vals):
+ print 'current prod: ', prod
+ # e.g. prod = ('20170505', '0', '12')
+ # ( date ,time, step)
+ # per date e.g. time = 0, 600, 1200, 1800
+ # per time e.g. step = 0, 3, 6, 9, 12
+ for i in range(len(index_keys)):
+ grib_index_select(iid, index_keys[i], prod[i])
+
+ gid = grib_new_from_index(iid)
+ # do convert2 program if gid at this time is not None,
+ # therefore save in hid
+ hid = gid
+ if gid is not None:
+ for k, f in fdict.iteritems():
+ silentremove(c.inputdir + "/fort." + k)
+ fdict[k] = open(c.inputdir + '/fort.' + k, 'w')
+
+ cdate = str(grib_get(gid, 'date'))
+ time = grib_get(gid, 'time')
+ type = grib_get(gid, 'type')
+ step = grib_get(gid, 'step')
+ # create correct timestamp from the three time informations
+ # date, time, step
+ timestamp = datetime.datetime.strptime(
+ cdate + '{:0>2}'.format(time/100), '%Y%m%d%H')
+ timestamp += datetime.timedelta(hours=int(step))
+
+ cdateH = datetime.datetime.strftime(timestamp, '%Y%m%d%H')
+ chms = datetime.datetime.strftime(timestamp, '%H%M%S')
+
+ if c.basetime is not None:
+ slimit = datetime.datetime.strptime(
+ c.start_date + '00', '%Y%m%d%H')
+ bt = '23'
+ if c.basetime == '00':
+ bt = '00'
+ slimit = datetime.datetime.strptime(
+ c.end_date + bt, '%Y%m%d%H') - \
+ datetime.timedelta(hours=12-int(c.dtime))
+ if c.basetime == '12':
+ bt = '12'
+ slimit = datetime.datetime.strptime(
+ c.end_date + bt, '%Y%m%d%H') - \
+ datetime.timedelta(hours=12-int(c.dtime))
+
+ elimit = datetime.datetime.strptime(
+ c.end_date + bt, '%Y%m%d%H')
+
+ if timestamp < slimit or timestamp > elimit:
+ continue
+
+ try:
+ if c.wrf == '1':
+ if 'olddate' not in locals():
+ fwrf = open(c.outputdir + '/WRF' + cdate +
+ '.{:0>2}'.format(time) + '.000.grb2', 'w')
+ olddate = cdate[:]
+ else:
+ if cdate != olddate:
+ fwrf = open(c.outputdir + '/WRF' + cdate +
+ '.{:0>2}'.format(time) + '.000.grb2',
+ 'w')
+ olddate = cdate[:]
+ except AttributeError:
+ pass
+
+
+ savedfields = []
+ while 1:
+ if gid is None:
+ break
+ paramId = grib_get(gid, 'paramId')
+ gridtype = grib_get(gid, 'gridType')
+ datatype = grib_get(gid, 'dataType')
+ levtype = grib_get(gid, 'typeOfLevel')
+ if paramId == 133 and gridtype == 'reduced_gg':
+ # Relative humidity (Q.grb) is used as a template only
+ # so we need the first we "meet"
+ with open(c.inputdir + '/fort.18', 'w') as fout:
+ grib_write(gid, fout)
+ elif paramId == 131 or paramId == 132:
+ grib_write(gid, fdict['10'])
+ elif paramId == 130:
+ grib_write(gid, fdict['11'])
+ elif paramId == 133 and gridtype != 'reduced_gg':
+ grib_write(gid, fdict['17'])
+ elif paramId == 152:
+ grib_write(gid, fdict['12'])
+ elif paramId == 155 and gridtype == 'sh':
+ grib_write(gid, fdict['13'])
+ elif paramId in [129, 138, 155] and levtype == 'hybrid' \
+ and c.wrf == '1':
+ pass
+ elif paramId == 246 or paramId == 247:
+ # cloud liquid water and ice
+ if paramId == 246:
+ clwc = grib_get_values(gid)
+ else:
+ clwc += grib_get_values(gid)
+ grib_set_values(gid, clwc)
+ grib_set(gid, 'paramId', 201031)
+ grib_write(gid, fdict['22'])
+ elif paramId == 135:
+ grib_write(gid, fdict['19'])
+ elif paramId == 77:
+ grib_write(gid, fdict['21'])
+ else:
+ if paramId not in savedfields:
+ grib_write(gid, fdict['16'])
+ savedfields.append(paramId)
+ else:
+ print('duplicate ' + str(paramId) + ' not written')
+
+ try:
+ if c.wrf == '1':
+# die if abfrage scheint ueberfluessig da eh das gleihce ausgefuehrt wird
+ if levtype == 'hybrid':
+ if paramId in [129, 130, 131, 132, 133, 138, 155]:
+ grib_write(gid, fwrf)
+ else:
+ if paramId in wrfpars:
+ grib_write(gid, fwrf)
+ except AttributeError:
+ pass
+
+ grib_release(gid)
+ gid = grib_new_from_index(iid)
+
+ for f in fdict.values():
+ f.close()
+
+ # call for CONVERT2
+# AUSLAGERN IN EIGENE FUNKTION
+
+ if hid is not None:
+ pwd = os.getcwd()
+ os.chdir(c.inputdir)
+ if os.stat('fort.21').st_size == 0 and int(c.eta) == 1:
+ print('Parameter 77 (etadot) is missing, most likely it is \
+ not available for this type or date/time\n')
+ print('Check parameters CLASS, TYPE, STREAM, START_DATE\n')
+ myerror(c, 'fort.21 is empty while parameter eta is set \
+ to 1 in CONTROL file')
+
+ p = subprocess.check_call([os.path.expandvars(
+ os.path.expanduser(c.exedir)) + '/CONVERT2'], shell=True)
+ os.chdir(pwd)
+ # create the corresponding output file fort.15
+ # (generated by CONVERT2)
+ # + fort.16 (paramId 167 and paramId 168)
+ fnout = c.inputdir + '/' + c.prefix
+ if c.maxstep > 12:
+ suffix = cdate[2:8] + '.{:0>2}'.format(time/100) + \
+ '.{:0>3}'.format(step)
+ else:
+ suffix = cdateH[2:10]
+
+ fnout += suffix
+ print("outputfile = " + fnout)
+ self.outputfilelist.append(fnout) # needed for final processing
+ fout = open(fnout, 'wb')
+ shutil.copyfileobj(open(c.inputdir + '/fort.15', 'rb'), fout)
+ if c.cwc == '1':
+ shutil.copyfileobj(open(c.inputdir + '/fort.22', 'rb'), fout)
+ shutil.copyfileobj(open(c.inputdir + '/flux' + cdate[0:2] +
+ suffix, 'rb'), fout)
+ shutil.copyfileobj(open(c.inputdir + '/fort.16', 'rb'), fout)
+ orolsm = glob.glob(c.inputdir +
+ '/OG_OROLSM__SL.*.' + c.ppid + '*')[0]
+ shutil.copyfileobj(open(orolsm, 'rb'), fout)
+ fout.close()
+ if c.omega == '1':
+ fnout = c.outputdir + '/OMEGA'
+ fout = open(fnout, 'wb')
+ shutil.copyfileobj(open(c.inputdir + '/fort.25', 'rb'), fout)
+
+ try:
+ if c.wrf == '1':
+ fwrf.close()
+ except:
+ pass
+
+ grib_index_release(iid)
+
+ return
+
+
+ def deacc_fluxes(self, inputfiles, c):
+ '''
+ @Description:
+
+
+ @Input:
+ self: instance of ECFlexpart
+ The current object of the class.
+
+ inputfiles: instance of UIOFiles
+ Contains a list of files.
+
+ c: instance of class Control
+ Contains all the parameters of control files, which are e.g.:
+ DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
+ STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
+ LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
+ OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
+ ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
+ MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
+ DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
+
+ For more information about format and content of the parameter
+ see documentation.
+
+ @Return:
+ <nothing>
+ '''
+
+ table128 = init128(c.ecmwfdatadir +
+ '/grib_templates/ecmwf_grib1_table_128')
+ pars = toparamId(self.params['OG_acc_SL'][0], table128)
+ index_keys = ["date", "time", "step"]
+ indexfile = c.inputdir + "/date_time_stepRange.idx"
+ silentremove(indexfile)
+ grib = GribTools(inputfiles.files)
+ # creates new index file
+ iid = grib.index(index_keys=index_keys, index_file=indexfile)
+
+ # read values of index keys
+ index_vals = []
+ for key in index_keys:
+ key_vals = grib_index_get(iid,key)
+ print(key_vals)
+ # have to sort the steps for disaggregation,
+ # therefore convert to int first
+ if key == 'step':
+ key_vals = [int(k) for k in key_vals]
+ key_vals.sort()
+ key_vals = [str(k) for k in key_vals]
+ index_vals.append(key_vals)
+ # index_vals looks for example like:
+ # index_vals[0]: ('20171106', '20171107', '20171108') ; date
+ # index_vals[1]: ('0', '1200') ; time
+ # index_vals[2]: (3', '6', '9', '12') ; stepRange
+
+ valsdict = {}
+ svalsdict = {}
+ stepsdict = {}
+ for p in pars:
+ valsdict[str(p)] = []
+ svalsdict[str(p)] = []
+ stepsdict[str(p)] = []
+
+ for prod in product(*index_vals):
+ # e.g. prod = ('20170505', '0', '12')
+ # ( date ,time, step)
+ # per date e.g. time = 0, 1200
+ # per time e.g. step = 3, 6, 9, 12
+ for i in range(len(index_keys)):
+ grib_index_select(iid, index_keys[i], prod[i])
+
+ gid = grib_new_from_index(iid)
+ # do convert2 program if gid at this time is not None,
+ # therefore save in hid
+ hid = gid
+ if gid is not None:
+ cdate = grib_get(gid, 'date')
+ time = grib_get(gid, 'time')
+ type = grib_get(gid, 'type')
+ step = grib_get(gid, 'step')
+ # date+time+step-2*dtime
+ # (since interpolated value valid for step-2*dtime)
+ sdate = datetime.datetime(year=cdate/10000,
+ month=mod(cdate, 10000)/100,
+ day=mod(cdate, 100),
+ hour=time/100)
+ fdate = sdate + datetime.timedelta(
+ hours=step-2*int(c.dtime))
+ sdates = sdate + datetime.timedelta(hours=step)
+ else:
+ break
+
+ if c.maxstep > 12:
+ fnout = c.inputdir + '/flux' + \
+ sdate.strftime('%Y%m%d') + '.{:0>2}'.format(time/100) + \
+ '.{:0>3}'.format(step-2*int(c.dtime))
+ gnout = c.inputdir + '/flux' + \
+ sdate.strftime('%Y%m%d') + '.{:0>2}'.format(time/100) + \
+ '.{:0>3}'.format(step-int(c.dtime))
+ hnout = c.inputdir + '/flux' + \
+ sdate.strftime('%Y%m%d') + '.{:0>2}'.format(time/100) + \
+ '.{:0>3}'.format(step)
+ g = open(gnout, 'w')
+ h = open(hnout, 'w')
+ else:
+ fnout = c.inputdir + '/flux' + fdate.strftime('%Y%m%d%H')
+ gnout = c.inputdir + '/flux' + (fdate+datetime.timedelta(
+ hours = int(c.dtime))).strftime('%Y%m%d%H')
+ hnout = c.inputdir + '/flux' + sdates.strftime('%Y%m%d%H')
+ g = open(gnout, 'w')
+ h = open(hnout, 'w')
+ print("outputfile = " + fnout)
+ f = open(fnout, 'w')
+
+ # read message for message and store relevant data fields
+ # data keywords are stored in pars
+ while 1:
+ if gid is None:
+ break
+ cparamId = str(grib_get(gid, 'paramId'))
+ step = grib_get(gid, 'step')
+ atime = grib_get(gid, 'time')
+ ni = grib_get(gid, 'Ni')
+ nj = grib_get(gid, 'Nj')
+ if cparamId in valsdict.keys():
+ values = grib_get_values(gid)
+ vdp = valsdict[cparamId]
+ svdp = svalsdict[cparamId]
+ sd = stepsdict[cparamId]
+
+ if cparamId == '142' or cparamId == '143':
+ fak = 1. / 1000.
+ else:
+ fak = 3600.
+
+ values = (reshape(values, (nj, ni))).flatten() / fak
+ vdp.append(values[:]) # save the accumulated values
+ if step <= int(c.dtime):
+ svdp.append(values[:] / int(c.dtime))
+ else: # deaccumulate values
+ svdp.append((vdp[-1] - vdp[-2]) / int(c.dtime))
+
+ print(cparamId, atime, step, len(values),
+ values[0], std(values))
+ # save the 1/3-hourly or specific values
+ # svdp.append(values[:])
+ sd.append(step)
+ # len(svdp) correspond to the time
+ if len(svdp) >= 3:
+ if len(svdp) > 3:
+ if cparamId == '142' or cparamId == '143':
+ values = Disagg.darain(svdp)
+ else:
+ values = Disagg.dapoly(svdp)
+
+ if not (step == c.maxstep and c.maxstep > 12 \
+ or sdates == elimit):
+ vdp.pop(0)
+ svdp.pop(0)
+ else:
+ if c.maxstep > 12:
+ values = svdp[1]
+ else:
+ values = svdp[0]
+
+ grib_set_values(gid, values)
+ if c.maxstep > 12:
+ grib_set(gid, 'step', max(0, step-2*int(c.dtime)))
+ else:
+ grib_set(gid, 'step', 0)
+ grib_set(gid, 'time', fdate.hour*100)
+ grib_set(gid, 'date', fdate.year*10000 +
+ fdate.month*100+fdate.day)
+ grib_write(gid, f)
+
+ if c.basetime is not None:
+ elimit = datetime.datetime.strptime(c.end_date +
+ c.basetime,
+ '%Y%m%d%H')
+ else:
+ elimit = sdate + datetime.timedelta(2*int(c.dtime))
+
+ # squeeze out information of last two steps contained
+ # in svdp
+ # if step+int(c.dtime) == c.maxstep and c.maxstep>12
+ # or sdates+datetime.timedelta(hours = int(c.dtime))
+ # >= elimit:
+ # Note that svdp[0] has not been popped in this case
+
+ if (step == c.maxstep and c.maxstep > 12
+ or sdates == elimit):
+
+ values = svdp[3]
+ grib_set_values(gid, values)
+ grib_set(gid, 'step', 0)
+ truedatetime = fdate + datetime.timedelta(
+ hours=2*int(c.dtime))
+ grib_set(gid, 'time', truedatetime.hour * 100)
+ grib_set(gid, 'date', truedatetime.year * 10000 +
+ truedatetime.month * 100 +
+ truedatetime.day)
+ grib_write(gid, h)
+
+ #values = (svdp[1]+svdp[2])/2.
+ if cparamId == '142' or cparamId == '143':
+ values = Disagg.darain(list(reversed(svdp)))
+ else:
+ values = Disagg.dapoly(list(reversed(svdp)))
+
+ grib_set(gid, 'step',0)
+ truedatetime = fdate + datetime.timedelta(
+ hours=int(c.dtime))
+ grib_set(gid, 'time', truedatetime.hour * 100)
+ grib_set(gid, 'date', truedatetime.year * 10000 +
+ truedatetime.month * 100 +
+ truedatetime.day)
+ grib_set_values(gid, values)
+ grib_write(gid, g)
+
+ grib_release(gid)
+
+ gid = grib_new_from_index(iid)
+
+ f.close()
+ g.close()
+ h.close()
+
+ grib_index_release(iid)
+
+ return
diff --git a/python/ECMWF_ENV b/python/ECMWF_ENV
index 5221e04..6706466 100644
--- a/python/ECMWF_ENV
+++ b/python/ECMWF_ENV
@@ -1,4 +1,4 @@
-ECUID lh0
-ECGID spatlh00
-GATEWAY srvx7.img.univie.ac.at
-DESTINATION leo@genericSftp
+ECUID km4a
+ECGID at
+GATEWAY srvx8.img.univie.ac.at
+DESTINATION philipa8@genericSftp
diff --git a/python/FlexpartTools.py b/python/FlexpartTools.py
deleted file mode 100644
index 27d3ad3..0000000
--- a/python/FlexpartTools.py
+++ /dev/null
@@ -1,2370 +0,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-#************************************************************************
-# TODO AP
-#AP
-# - Functionality Description Provided is not correct for this file
-# - localpythonpath should not be set in module load section!
-# - Change History ist nicht angepasst ans File! Original geben lassen
-# - def myerror muss angepasst werden da derzeit manuelle modifikation notwendig
-# - the gaussian keyword in mars retrieval should be removed!!! but kept in
-# control file for the reason of calculation of eta dot from gaussian grid
-# - call of convert in eigene Funktion auslagern
-
-#************************************************************************
-"""
-@Author: Anne Fouilloux (University of Oslo)
-
-@Date: October 2014
-
-@ChangeHistory:
- November 2015 - Leopold Haimberger (University of Vienna):
- - using the WebAPI also for general MARS retrievals
- - job submission on ecgate and cca
- - job templates suitable for twice daily operational dissemination
- - dividing retrievals of longer periods into digestable chunks
- - retrieve also longer term forecasts, not only analyses and
- short term forecast data
- - conversion into GRIB2
- - conversion into .fp format for faster execution of FLEXPART
-
- February 2018 - Anne Philipp (University of Vienna):
- - applied PEP8 style guide
- - added documentation
-
-@License:
- (C) Copyright 2014 UIO.
-
- This software is licensed under the terms of the Apache Licence Version 2.0
- which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
-
-@Requirements:
- - A standard python 2.6 or 2.7 installation
- - dateutils
- - matplotlib (optional, for debugging)
- - ECMWF specific packages, all available from https://software.ecmwf.int/
- ECMWF WebMARS, gribAPI with python enabled, emoslib and
- ecaccess web toolkit
-
-@Description:
- Further documentation may be obtained from www.flexpart.eu.
-
- Functionality provided:
- Prepare input 3D-wind fields in hybrid coordinates +
- surface fields for FLEXPART runs
-"""
-# ------------------------------------------------------------------------------
-# MODULES
-# ------------------------------------------------------------------------------
-import subprocess
-from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter
-import traceback
-import shutil
-import os
-import errno
-import sys
-import inspect
-import glob
-import datetime
-from string import atoi
-from numpy import *
-ecapi = True
-try:
- import ecmwfapi
-except ImportError:
- ecapi = False
-from gribapi import *
-from GribTools import GribTools
-
-def interpret_args_and_control():
- '''
- @Description:
- Assigns the command line arguments and reads control file
- content. Apply default values for non mentioned arguments.
-
- @Input:
- <nothing>
-
- @Return:
- args: instance of ArgumentParser
- Contains the commandline arguments from script/program call.
-
- c: instance of class Control
- Contains all necessary information of a control file. The parameters
- are: DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM,
- NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL, LEVELIST,
- RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA,
- SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS,
- ECFSDIR, MAILOPS, MAILFAIL, GRIB2FLEXPART, DEBUG, INPUTDIR,
- OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
- For more information about format and content of the parameter see
- documentation.
-
- '''
- parser = ArgumentParser(description='Retrieve FLEXPART input from \
- ECMWF MARS archive',
- formatter_class=ArgumentDefaultsHelpFormatter)
-
- # the most important arguments
- parser.add_argument("--start_date", dest="start_date",
- help="start date YYYYMMDD")
- parser.add_argument("--end_date", dest="end_date",
- help="end_date YYYYMMDD")
- parser.add_argument("--date_chunk", dest="date_chunk", default=None,
- help="# of days to be retrieved at once")
-
- # some arguments that override the default in the control file
- parser.add_argument("--basetime", dest="basetime",
- help="base such as 00/12 (for half day retrievals)")
- parser.add_argument("--step", dest="step",
- help="steps such as 00/to/48")
- parser.add_argument("--levelist", dest="levelist",
- help="Vertical levels to be retrieved, e.g. 30/to/60")
- parser.add_argument("--area", dest="area",
- help="area defined as north/west/south/east")
-
- # set the working directories
- parser.add_argument("--inputdir", dest="inputdir", default=None,
- help="root directory for storing intermediate files")
- parser.add_argument("--outputdir", dest="outputdir", default=None,
- help="root directory for storing output files")
- parser.add_argument("--flexpart_root_scripts", dest="flexpart_root_scripts",
- help="FLEXPART root directory (to find grib2flexpart \
- and COMMAND file)\n\ Normally ECMWFDATA resides in \
- the scripts directory of the FLEXPART distribution")
-
- # this is only used by prepareFLEXPART.py to rerun a postprocessing step
- parser.add_argument("--ppid", dest="ppid",
- help="Specify parent process id for \
- rerun of prepareFLEXPART")
-
- # arguments for job submission to ECMWF, only needed by submit.py
- parser.add_argument("--job_template", dest='job_template',
- default="job.temp",
- help="job template file for submission to ECMWF")
- parser.add_argument("--queue", dest="queue",
- help="queue for submission to ECMWF \
- (e.g. ecgate or cca )")
- parser.add_argument("--controlfile", dest="controlfile",
- default='CONTROL.temp',
- help="file with control parameters")
- parser.add_argument("--debug", dest="debug", default=0,
- help="Debug mode - leave temporary files intact")
-
- args = parser.parse_args()
-
- # create instance of Control for specified controlfile
- # and assign the parameters (and default values if necessary)
- try:
- c = Control(args.controlfile)
- except IOError:
- try:
- c = Control(localpythonpath + args.controlfile)
- except:
- print('Could not read control file "' + args.controlfile + '"')
- print('Either it does not exist or its syntax is wrong.')
- print('Try "' + sys.argv[0].split('/')[-1] +
- ' -h" to print usage information')
- exit(1)
-
- # check for having at least a starting date
- if args.start_date is None and getattr(c, 'start_date') is None:
- print('start_date specified neither in command line nor \
- in control file ' + args.controlfile)
- print('Try "' + sys.argv[0].split('/')[-1] +
- ' -h" to print usage information')
- exit(1)
-
- # save all existing command line parameter to the Control instance
- # if parameter is not specified through the command line or CONTROL file
- # set default values
- if args.start_date is not None:
- c.start_date = args.start_date
- if args.end_date is not None:
- c.end_date = args.end_date
- if c.end_date is None:
- c.end_date = c.start_date
- if args.date_chunk is not None:
- c.date_chunk = args.date_chunk
-
- if not hasattr(c, 'debug'):
- c.debug = args.debug
-
- if args.inputdir is None and args.outputdir is None:
- c.inputdir = '../work'
- c.outputdir = '../work'
- else:
- if args.inputdir is not None:
- c.inputdir = args.inputdir
- if args.outputdir is None:
- c.outputdir = args.inputdir
- if args.outputdir is not None:
- c.outputdir = args.outputdir
- if args.inputdir is None:
- c.inputdir = args.outputdir
-
- if hasattr(c, 'outputdir') is False and args.outputdir is None:
- c.outputdir = c.inputdir
- else:
- if args.outputdir is not None:
- c.outputdir = args.outputdir
-
- if args.area is not None:
- afloat = '.' in args.area
- l = args.area.split('/')
- if afloat:
- for i in range(len(l)):
- l[i] = str(int(float(l[i]) * 1000))
- c.upper, c.left, c.lower, c.right = l
-
- # NOTE: basetime activates the ''operational mode''
- if args.basetime is not None:
- c.basetime = args.basetime
-
- if args.step is not None:
- l = args.step.split('/')
- if 'to' in args.step.lower():
- if 'by' in args.step.lower():
- ilist = arange(int(l[0]), int(l[2]) + 1, int(l[4]))
- c.step = ['{:0>3}'.format(i) for i in ilist]
- else:
- myerror(None, args.step + ':\n' +
- 'please use "by" as well if "to" is used')
- else:
- c.step = l
-
- if args.levelist is not None:
- c.levelist = args.levelist
- if 'to' in c.levelist:
- c.level = c.levelist.split('/')[2]
- else:
- c.level = c.levelist.split('/')[-1]
-
- if args.flexpart_root_scripts is not None:
- c.flexpart_root_scripts = args.flexpart_root_scripts
-
- return args, c
-
-
-def install_args_and_control():
- '''
- @Description:
- Assigns the command line arguments for installation and reads
- control file content. Apply default values for non mentioned arguments.
-
- @Input:
- <nothing>
-
- @Return:
- args: instance of ArgumentParser
- Contains the commandline arguments from script/program call.
-
- c: instance of class Control
- Contains all necessary information of a control file. The parameters
- are: DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM,
- NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL, LEVELIST,
- RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA,
- SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS,
- ECFSDIR, MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR
- For more information about format and content of the parameter see
- documentation.
- '''
- parser = ArgumentParser(description='Install ECMWFDATA software locally or \
- on ECMWF machines',
- formatter_class=ArgumentDefaultsHelpFormatter)
-
- parser.add_argument('--target', dest='install_target',
- help="Valid targets: local | ecgate | cca , \
- the latter two are at ECMWF")
- parser.add_argument("--makefile", dest="makefile",
- help='Name of Makefile to use for compiling CONVERT2')
- parser.add_argument("--ecuid", dest="ecuid",
- help='user id at ECMWF')
- parser.add_argument("--ecgid", dest="ecgid",
- help='group id at ECMWF')
- parser.add_argument("--gateway", dest="gateway",
- help='name of local gateway server')
- parser.add_argument("--destination", dest="destination",
- help='ecaccess destination, e.g. leo@genericSftp')
-
- parser.add_argument("--flexpart_root_scripts", dest="flexpart_root_scripts",
- help="FLEXPART root directory on ECMWF servers \
- (to find grib2flexpart and COMMAND file)\n\
- Normally ECMWFDATA resides in the scripts directory \
- of the FLEXPART distribution, thus the:")
-
-# arguments for job submission to ECMWF, only needed by submit.py
- parser.add_argument("--job_template", dest='job_template',
- default="job.temp.o",
- help="job template file for submission to ECMWF")
-
- parser.add_argument("--controlfile", dest="controlfile",
- default='CONTROL.temp',
- help="file with control parameters")
-
- args = parser.parse_args()
-
- try:
- c = Control(args.controlfile)
- except:
- print('Could not read control file "' + args.controlfile + '"')
- print('Either it does not exist or its syntax is wrong.')
- print('Try "' + sys.argv[0].split('/')[-1] +
- ' -h" to print usage information')
- exit(1)
-
- if args.install_target != 'local':
- if (args.ecgid is None or args.ecuid is None or args.gateway is None
- or args.destination is None):
- print('Please enter your ECMWF user id and group id as well as \
- the \nname of the local gateway and the ectrans \
- destination ')
- print('with command line options --ecuid --ecgid \
- --gateway --destination')
- print('Try "' + sys.argv[0].split('/')[-1] +
- ' -h" to print usage information')
- print('Please consult ecaccess documentation or ECMWF user support \
- for further details')
- sys.exit(1)
- else:
- c.ecuid = args.ecuid
- c.ecgid = args.ecgid
- c.gateway = args.gateway
- c.destination = args.destination
-
- try:
- c.makefile = args.makefile
- except:
- pass
-
- if args.install_target == 'local':
- if args.flexpart_root_scripts is None:
- c.flexpart_root_scripts = '../'
- else:
- c.flexpart_root_scripts = args.flexpart_root_scripts
-
- if args.install_target != 'local':
- if args.flexpart_root_scripts is None:
- c.ec_flexpart_root_scripts = '${HOME}'
- else:
- c.ec_flexpart_root_scripts = args.flexpart_root_scripts
-
- return args, c
-
-
-def cleanup(c):
- '''
- @Description:
- Remove all files from intermediate directory
- (inputdir from control file).
-
- @Input:
- c: instance of class Control
- Contains all the parameters of control files, which are e.g.:
- DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
- STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
- LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
- OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
- ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
- MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
- DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
-
- For more information about format and content of the parameter
- see documentation.
-
- @Return:
- <nothing>
- '''
-
- print("cleanup")
-
- cleanlist = glob.glob(c.inputdir + "/*")
- for cl in cleanlist:
- if c.prefix not in cl:
- silentremove(cl)
- if c.ecapi is False and (c.ectrans == '1' or c.ecstorage == '1'):
- silentremove(cl)
-
- print("Done")
-
- return
-
-
-def myerror(c, message='ERROR'):
- '''
- @Description:
- Prints a specified error message which can be passed to the function
- before exiting the program.
-
- @Input:
- c: instance of class Control
- Contains all the parameters of control files, which are e.g.:
- DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
- STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
- LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
- OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
- ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
- MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
- DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
-
- For more information about format and content of the parameter
- see documentation.
-
- message: string, optional
- Error message. Default value is "ERROR".
-
- @Return:
- <nothing>
- '''
- # uncomment if user wants email notification directly from python
- #try:
- #target = c.mailfail
- #except AttributeError:
- #target = os.getenv('USER')
-
- #if(type(target) is not list):
- #target = [target]
-
- print(message)
-
- # uncomment if user wants email notification directly from python
- #for t in target:
- #p = subprocess.Popen(['mail','-s ECMWFDATA v7.0 ERROR', os.path.expandvars(t)],
- # stdin = subprocess.PIPE, stdout = subprocess.PIPE,
- # stderr = subprocess.PIPE, bufsize = 1)
- #tr = '\n'.join(traceback.format_stack())
- #pout = p.communicate(input = message+'\n\n'+tr)[0]
- #print 'Email sent to '+os.path.expandvars(t) # +' '+pout.decode()
-
- exit(1)
-
- return
-
-
-def normalexit(c, message='Done!'):
- '''
- @Description:
- Prints a specific exit message which can be passed to the function.
-
- @Input:
- c: instance of class Control
- Contains all the parameters of control files, which are e.g.:
- DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
- STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
- LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
- OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
- ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
- MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
- DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
-
- For more information about format and content of the parameter
- see documentation.
-
- message: string, optional
- Message for exiting program. Default value is "Done!".
-
- @Return:
- <nothing>
-
- '''
- # Uncomment if user wants notification directly from python
- #try:
- #target = c.mailops
- #if(type(target) is not list):
- #target = [target]
- #for t in target:
- #p = subprocess.Popen(['mail','-s ECMWFDATA v7.0 normal exit',
- # os.path.expandvars(t)],
- # stdin = subprocess.PIPE,
- # stdout = subprocess.PIPE,
- # stderr = subprocess.PIPE, bufsize = 1)
- #pout = p.communicate(input = message+'\n\n')[0]
- #print pout.decode()
- #except:
- #pass
-
- print(message)
-
- return
-
-
-def product(*args, **kwds):
- '''
- @Description:
- This method is taken from an example at the ECMWF wiki website.
- https://software.ecmwf.int/wiki/display/GRIB/index.py; 2018-03-16
-
- This method combines the single characters of the passed arguments
- with each other. So that each character of each argument value
- will be combined with each character of the other arguments as a tuple.
-
- Example:
- product('ABCD', 'xy') --> Ax Ay Bx By Cx Cy Dx Dy
- product(range(2), repeat = 3) --> 000 001 010 011 100 101 110 111
-
- @Input:
- *args: tuple
- Positional arguments (arbitrary number).
-
- **kwds: dictionary
- Contains all the keyword arguments from *args.
-
- @Return:
- prod: tuple
- Return will be done with "yield". A tuple of combined arguments.
- See example in description above.
- '''
-
- pools = map(tuple, args) * kwds.get('repeat', 1)
- result = [[]]
- for pool in pools:
- result = [x + [y] for x in result for y in pool]
- for prod in result:
- yield tuple(prod)
-
- return
-
-
-def silentremove(filename):
- '''
- @Description:
- Removes the file which name is passed to the function if
- it exists. The function does not fail if the file does not
- exist.
-
- @Input:
- filename: string
- The name of the file to be removed without notification.
-
- @Return:
- <nothing>
- '''
- try:
- os.remove(filename)
- except OSError as e:
- # this would be "except OSError, e:" before Python 2.6
- if e.errno is not errno.ENOENT:
- # errno.ENOENT = no such file or directory
- raise # re-raise exception if a different error occured
-
- return
-
-
-def init128(fn):
- '''
- @Description:
- Opens and reads the grib file with table 128 information.
-
- @Input:
- fn: string
- Path to file of ECMWF grib table number 128.
-
- @Return:
- table128: dictionary
- Contains the ECMWF grib table 128 information.
- '''
- table128 = dict()
- with open(fn) as f:
- fdata = f.read().split('\n')
- for data in fdata:
- if data[0] != '!':
- table128[data[0:3]] = data[59:64].strip()
-
- return table128
-
-
-def toparamId(pars, table):
- '''
- @Description:
- Transform parameter names to parameter ids
- with ECMWF grib table 128.
-
- @Input:
- pars: string
- Addpar argument from control file in the format of
- parameter names instead of ids.
-
- table: dictionary
- Contains the ECMWF grib table 128 information.
-
- @Return:
- ipar: list of integer
- List of addpar parameters from control file transformed to
- parameter ids in the format of integer.
- '''
- cpar = pars.upper().split('/')
- ipar = []
- for par in cpar:
- found = False
- for k, v in table.iteritems():
- if par == k or par == v:
- ipar.append(int(k))
- found = True
- break
- if found is False:
- print('Warning: par ' + par + ' not found in table 128')
-
- return ipar
-
-
-def dapoly(alist):
- '''
- @Description:
- Interpolation of deaccumulated fluxes of an ECMWF model FG field
- using a cubic polynomial solution which conserves the integrals
- of the fluxes within each timespan.
- Disaggregation is done for a list of 4 values to generate a new
- interpolated value which is output at the central point of the 4
- accumulation timespans.
-
- @Input:
- alist: list of size 4, float
- List of 4 flux values.
-
- @Return:
- nvalue: float
- New value for the second position of the disaggregated
- fluxes field.
-
- '''
- pya = (alist[3] - alist[0] + 3. * (alist[1] - alist[2])) / 6.
- pyb = (alist[2] + alist[0]) / 2. - alist[1] - 9. * pya / 2.
- pyc = alist[1] - alist[0] - 7. * pya / 2. - 2. * pyb
- pyd = alist[0] - pya / 4. - pyb / 3. - pyc / 2.
- nvalue = 8. * pya + 4. * pyb + 2. * pyc + pyd
-
- return nvalue
-
-
-def darain(alist):
- '''
- @Description:
- Interpolation of deaccumulated precipitation fiels of the ECMWF fields
- using a modified linear solution.
- Disaggregate a list of 4 precipitation values to generate a new value
- for the second position of the 4 value list. The interpolated values
- are output at the central point of the 4 accumulation timespans
- This is used for precipitation fields.
-
- @Input:
- alist: list of size 4, float
- List of 4 precipitation values.
-
- @Return:
- nvalue: float
- New value for the second position of the disaggregated
- precipitation field.
- '''
- xa = alist[0]
- xb = alist[1]
- xc = alist[2]
- xd = alist[3]
- xa[xa < 0.] = 0.
- xb[xb < 0.] = 0.
- xc[xc < 0.] = 0.
- xd[xd < 0.] = 0.
-
- xac = 0.5 * xb
- mask = xa + xc > 0.
- xac[mask] = xb[mask] * xc[mask] / (xa[mask] + xc[mask])
- xbd = 0.5 * xc
- mask = xb + xd > 0.
- xbd[mask] = xb[mask] * xc[mask] / (xb[mask] + xd[mask])
- nvalue = xac + xbd
-
- return nvalue
-
-
-class Control:
- '''
- Class containing the information of the ECMWFDATA control file.
-
- Contains all the parameters of control files, which are e.g.:
- DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
- STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
- LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
- OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
- ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
- MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR,
- BASETIME, DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
-
- For more information about format and content of the parameter
- see documentation.
-
- '''
-
- def __init__(self, filename):
- '''
- @Description:
- Initialises the instance of Control class and defines and
- assign all controlfile variables. Set default values if
- parameter was not in CONTROL file.
-
- @Input:
- self: instance of Control class
- Description see class documentation.
-
- filename: string
- Name of control file.
-
- @Return:
- <nothing>
- '''
-
- # read whole CONTROL file
- with open(filename) as f:
- fdata = f.read().split('\n')
-
- # go through every line and store parameter
- # as class variable
- for ldata in fdata:
- data = ldata.split()
- if len(data) > 1:
- if 'm_' in data[0].lower():
- data[0] = data[0][2:]
- if data[0].lower() == 'class':
- data[0] = 'marsclass'
- if data[0].lower() == 'day1':
- data[0] = 'start_date'
- if data[0].lower() == 'day2':
- data[0] = 'end_date'
- if data[0].lower() == 'addpar':
- if '/' in data[1]:
- # remove leading '/' sign from addpar content
- if data[1][0] == '/':
- data[1] = data[1][1:]
- dd = data[1].split('/')
- data = [data[0]]
- for d in dd:
- data.append(d)
- pass
- if len(data) == 2:
- if '$' in data[1]:
- setattr(self, data[0].lower(), data[1])
- while '$' in data[1]:
- i = data[1].index('$')
- j = data[1].find('{')
- k = data[1].find('}')
- var = os.getenv(data[1][j+1:k])
- if var is not None:
- data[1] = data[1][:i] + var + data[1][k+1:]
- else:
- myerror(None, 'Could not find variable ' +
- data[1][j+1:k] + ' while reading ' +
- filename)
- setattr(self, data[0].lower()+'_expanded', data[1])
- else:
- if data[1].lower() != 'none':
- setattr(self, data[0].lower(), data[1])
- else:
- setattr(self, data[0].lower(), None)
- elif len(data) > 2:
- setattr(self, data[0].lower(), (data[1:]))
- else:
- pass
-
- # check a couple of necessary attributes if they contain values
- # otherwise set default values
- if not hasattr(self, 'start_date'):
- self.start_date = None
- if not hasattr(self, 'end_date'):
- self.end_date = self.start_date
- if not hasattr(self, 'accuracy'):
- self.accuracy = 24
- if not hasattr(self, 'omega'):
- self.omega = '0'
- if not hasattr(self, 'cwc'):
- self.cwc = '0'
- if not hasattr(self, 'omegadiff'):
- self.omegadiff = '0'
- if not hasattr(self, 'etadiff'):
- self.etadiff = '0'
- if not hasattr(self, 'levelist'):
- if not hasattr(self, 'level'):
- print(('Warning: neither levelist nor level \
- specified in CONTROL file'))
- else:
- self.levelist = '1/to/' + self.level
- else:
- if 'to' in self.levelist:
- self.level = self.levelist.split('/')[2]
- else:
- self.level = self.levelist.split('/')[-1]
-
- if not hasattr(self, 'maxstep'):
- # find out maximum step
- self.maxstep = 0
- for s in self.step:
- if int(s) > self.maxstep:
- self.maxstep = int(s)
- else:
- self.maxstep = int(self.maxstep)
-
- if not hasattr(self, 'prefix'):
- self.prefix = 'EN'
- if not hasattr(self, 'makefile'):
- self.makefile = None
- if not hasattr(self, 'basetime'):
- self.basetime = None
- if not hasattr(self, 'date_chunk'):
- self.date_chunk = '3'
- if not hasattr(self, 'grib2flexpart'):
- self.grib2flexpart = '0'
-
- # script directory
- self.ecmwfdatadir = os.path.dirname(os.path.abspath(
- inspect.getfile(
- inspect.currentframe()))) + '/../'
- # Fortran source directory
- self.exedir = self.ecmwfdatadir + 'src/'
-
- # FLEXPART directory
- if not hasattr(self, 'flexpart_root_scripts'):
- self.flexpart_root_scripts = self.ecmwfdatadir
-
- return
-
-
- def __str__(self):
- '''
- @Description:
- Prepares a single string with all the comma seperated Control
- class attributes including their values.
-
- Example:
- {'kids': 0, 'name': 'Dog', 'color': 'Spotted',
- 'age': 10, 'legs': 2, 'smell': 'Alot'}
-
- @Input:
- self: instance of Control class
- Description see class documentation.
-
- @Return:
- string of Control class attributes with their values
- '''
-
- attrs = vars(self)
-
- return ', '.join("%s: %s" % item for item in attrs.items())
-
- def tolist(self):
- '''
- @Description:
- Just generates a list of strings containing the attributes and
- assigned values except the attributes "_expanded", "exedir",
- "ecmwfdatadir" and "flexpart_root_scripts".
-
- @Input:
- self: instance of Control class
- Description see class documentation.
-
- @Return:
- l: list
- A sorted list of the all Control class attributes with
- their values except the attributes "_expanded", "exedir",
- "ecmwfdatadir" and "flexpart_root_scripts".
- '''
- attrs = vars(self)
- l = list()
- for item in attrs.items():
- if '_expanded' in item[0]:
- pass
- elif 'exedir' in item[0]:
- pass
- elif 'flexpart_root_scripts' in item[0]:
- pass
- elif 'ecmwfdatadir' in item[0]:
- pass
- else:
- if type(item[1]) is list:
- stot = ''
- for s in item[1]:
- stot += s + ' '
-
- l.append("%s %s" % (item[0], stot))
- else:
-#AP syntax error with doubled %s ???
- l.append("%s %s" % item )
- return sorted(l)
-
-
-class MARSretrieval:
- '''
- Class for submitting MARS retrievals.
-
- A description of MARS keywords/arguments and examples of their
- values can be found here:
- https://software.ecmwf.int/wiki/display/UDOC/\
- Identification+keywords#Identificationkeywords-class
-
- '''
-
- def __init__(self, server, marsclass = "ei", type = "", levtype = "",
- levelist = "", repres = "", date = "", resol = "", stream = "",
- area = "", time = "", step = "", expver = "1", number = "",
- accuracy = "", grid = "", gaussian = "", target = "",
- param = ""):
- '''
- @Description:
- Initialises the instance of the MARSretrieval class and
- defines and assigns a set of the necessary retrieval parameters
- for the FLEXPART input data.
- A description of MARS keywords/arguments, their dependencies
- on each other and examples of their values can be found here:
-
- https://software.ecmwf.int/wiki/display/UDOC/MARS+keywords
-
- @Input:
- self: instance of MARSretrieval
- For description see class documentation.
-
- server: instance of ECMWFService (from ECMWF Web-API)
- This is the connection to the ECMWF data servers.
- It is needed for the pythonic access of ECMWF data.
-
- marsclass: string, optional
- Characterisation of dataset. E.g. EI (ERA-Interim),
- E4 (ERA40), OD (Operational archive), ea (ERA5).
- Default is the ERA-Interim dataset "ei".
-
- type: string, optional
- Determines the type of fields to be retrieved.
- Selects between observations, images or fields.
- Examples for fields: Analysis (an), Forecast (fc),
- Perturbed Forecast (pf), Control Forecast (cf) and so on.
- Default is an empty string.
-
- levtype: string, optional
- Denotes type of level. Has a direct implication on valid
- levelist values!
- E.g. model level (ml), pressure level (pl), surface (sfc),
- potential vorticity (pv), potential temperature (pt)
- and depth (dp).
- Default is an empty string.
-
- levelist: string, optional
- Specifies the required levels. It has to have a valid
- correspondence to the selected levtype.
- Examples: model level: 1/to/137, pressure levels: 500/to/1000
- Default is an empty string.
-
- repres: string, optional
- Selects the representation of the archived data.
- E.g. sh - spherical harmonics, gg - Gaussian grid,
- ll - latitude/longitude, ...
- Default is an empty string.
-
- date: string, optional
- Specifies the Analysis date, the Forecast base date or
- Observations date. Valid formats are:
- Absolute as YYYY-MM-DD or YYYYMMDD.
- Default is an empty string.
-
- resol: string, optional
- Specifies the desired triangular truncation of retrieved data,
- before carrying out any other selected post-processing.
- The default is automatic truncation (auto), by which the lowest
- resolution compatible with the value specified in grid is
- automatically selected for the retrieval.
- Users wanting to perform post-processing from full spectral
- resolution should specify Archived Value (av).
- The following are examples of existing resolutions found in
- the archive: 63, 106, 159, 213, 255, 319, 399, 511, 799 or 1279.
- This keyword has no meaning/effect if the archived data is
- not in spherical harmonics representation.
- The best selection can be found here:
- https://software.ecmwf.int/wiki/display/UDOC/\
- Retrieve#Retrieve-Truncationbeforeinterpolation
- Default is an empty string.
-
- stream: string, optional
- Identifies the forecasting system used to generate the data.
- E.g. oper (Atmospheric model), enfo (Ensemble forecats), ...
- Default is an empty string.
-
- area: string, optional
- Specifies the desired sub-area of data to be extracted.
- Areas can be defined to wrap around the globe.
-
- Latitude values must be given as signed numbers, with:
- north latitudes (i.e. north of the equator)
- being positive (e.g: 40.5)
- south latitutes (i.e. south of the equator)
- being negative (e.g: -50.5)
- Longtitude values must be given as signed numbers, with:
- east longitudes (i.e. east of the 0 degree meridian)
- being positive (e.g: 35.0)
- west longitudes (i.e. west of the 0 degree meridian)
- being negative (e.g: -20.5)
-
- E.g.: North/West/South/East
- Default is an empty string.
-
- time: string, optional
- Specifies the time of the data in hours and minutes.
- Valid values depend on the type of data: Analysis time,
- Forecast base time or First guess verification time
- (all usually at synoptic hours: 00, 06, 12 and 18 ).
- Observation time (any combination in hours and minutes is valid,
- subject to data availability in the archive).
- The syntax is HHMM or HH:MM. If MM is omitted it defaults to 00.
- Default is an empty string.
-
- step: string, optional
- Specifies the forecast time step from forecast base time.
- Valid values are hours (HH) from forecast base time. It also
- specifies the length of the forecast which verifies at
- First Guess time.
- E.g. 1/3/6-hourly
- Default is an empty string.
-
- expver: string, optional
- The version of the dataset. Each experiment is assigned a
- unique code (version). Production data is assigned 1 or 2,
- and experimental data in Operations 11, 12 ,...
- Research or Member State's experiments have a four letter
- experiment identifier.
- Default is "1".
-
- number: string, optional
- Selects the member in ensemble forecast run. (Only then it
- is necessary.) It has a different meaning depending on
- the type of data.
- E.g. Perturbed Forecasts: specifies the Ensemble forecast member
- Default is an empty string.
-
- accuracy: string, optional
- Specifies the number of bits per value to be used in the
- generated GRIB coded fields.
- A positive integer may be given to specify the preferred number
- of bits per packed value. This must not be greater than the
- number of bits normally used for a Fortran integer on the
- processor handling the request (typically 32 or 64 bit).
- Within a compute request the accuracy of the original fields
- can be passed to the result field by specifying accuracy=av.
- Default is an empty string.
-
- grid: string, optional
- Specifies the output grid which can be either a Gaussian grid
- or a Latitude/Longitude grid. MARS requests specifying
- grid=av will return the archived model grid.
-
- Lat/Lon grid: The grid spacing needs to be an integer
- fraction of 90 degrees e.g. grid = 0.5/0.5
-
- Gaussian grid: specified by a letter denoting the type of
- Gaussian grid followed by an integer (the grid number)
- representing the number of lines between the Pole and Equator,
- e.g.
- grid = F160 - full (or regular) Gaussian grid with
- 160 latitude lines between the pole and equator
- grid = N320 - ECMWF original reduced Gaussian grid with
- 320 latitude lines between the pole and equator,
- see Reduced Gaussian Grids for grid numbers used at ECMWF
- grid = O640 - ECMWF octahedral (reduced) Gaussian grid with
- 640 latitude lines between the pole and equator
- Default is an empty string.
-
- gaussian: string, optional
- This parameter is deprecated and should no longer be used.
- Specifies the desired type of Gaussian grid for the output.
- Valid Gaussian grids are quasi-regular (reduced) or regular.
- Keyword gaussian can only be specified together with
- keyword grid. Gaussian without grid has no effect.
- Default is an empty string.
-
- target: string, optional
- Specifies a file into which data is to be written after
- retrieval or manipulation. Path names should always be
- enclosed in double quotes. The MARS client supports automatic
- generation of multiple target files using MARS keywords
- enclosed in square brackets [ ]. If the environment variable
- MARS_MULTITARGET_STRICT_FORMAT is set to 1 before calling mars,
- the keyword values will be used in the filename as shown by
- the ecCodes GRIB tool grib_ls -m, e.g. with
- MARS_MULTITARGET_STRICT_FORMAT set to 1 the keywords time,
- expver and param will be formatted as 0600, 0001 and 129.128
- rather than 600, 1 and 129.
- Default is an empty string.
-
- param: string, optional
- Specifies the meteorological parameter.
- The list of meteorological parameters in MARS is extensive.
- Their availability is directly related to their meteorological
- meaning and, therefore, the rest of directives specified
- in the MARS request.
- Meteorological parameters can be specified by their
- GRIB code (param=130), their mnemonic (param=t) or
- full name (param=temperature).
- The list of parameter should be seperated by a "/"-sign.
- E.g. 130/131/133
- Default is an empty string.
-
- @Return:
- <nothing>
- '''
-
- self.server = server
- self.marsclass = marsclass
- self.type = type
- self.levtype = levtype
- self.levelist = levelist
- self.repres = repres
- self.date = date
- self.resol = resol
- self.stream = stream
- self.area = area
- self.time = time
- self.step = step
- self.expver = expver
- self.number = number
- self.accuracy = accuracy
- self.grid = grid
- self.gaussian = gaussian
- self.target = target
- self.param = param
-
- return
-
-
- def displayInfo(self):
- '''
- @Description:
- Prints all class attributes and their values.
-
- @Input:
- self: instance of MARSretrieval
- For description see class documentation.
-
- @Return:
- <nothing>
- '''
- # Get all class attributes and their values as a dictionary
- attrs = vars(self)
-
- # iterate through all attributes and print them
- # with their corresponding values
- for item in attrs.items():
- if item[0] in ('server'):
- pass
- else:
- print(item[0] + ': ' + str(item[1]))
-
- return
-
- def dataRetrieve(self):
- '''
- @Description:
- Submits a MARS retrieval. Depending on the existence of
- ECMWF Web-API it is submitted via Python or a
- subprocess in the Shell. The parameter for the mars retrieval
- are taken from the defined class attributes.
-
- @Input:
- self: instance of MARSretrieval
- For description see class documentation.
-
- @Return:
- <nothing>
- '''
- # Get all class attributes and their values as a dictionary
- attrs = vars(self)
-
- # convert the dictionary of attributes into a comma
- # seperated list of attributes with their values
- # needed for the retrieval call
- s = 'ret'
- for k, v in attrs.iteritems():
- if k in ('server'):
- continue
- if k == 'marsclass':
- k = 'class'
- if v == '':
- continue
- if k.lower() == 'target':
- target = v
- else:
- s = s + ',' + k + '=' + str(v)
-
- # MARS request via Python script
- if self.server is not False:
- try:
- self.server.execute(s, target)
- except:
- print('MARS Request failed, \
- have you already registered at apps.ecmwf.int?')
- raise IOError
- if os.stat(target).st_size == 0:
- print('MARS Request returned no data - please check request')
- raise IOError
- # MARS request via extra process in shell
- else:
- s += ',target = "' + target + '"'
- p = subprocess.Popen(['mars'], stdin=subprocess.PIPE,
- stdout=subprocess.PIPE,
- stderr=subprocess.PIPE, bufsize=1)
- pout = p.communicate(input=s)[0]
- print(pout.decode())
- if 'Some errors reported' in pout.decode():
- print('MARS Request failed - please check request')
- raise IOError
-
- if os.stat(target).st_size == 0:
- print('MARS Request returned no data - please check request')
- raise IOError
-
- return
-
-
-
-
-class ECFlexpart:
- '''
- Class to retrieve ECMWF data specific for running FLEXPART.
- '''
- def __init__(self, c, fluxes=False): #done/ verstehen
- '''
- @Description:
- Creates an object/instance of ECFlexpart with the
- associated settings of its attributes for the retrieval.
-
- @Input:
- self: instance of ECFlexpart
- The current object of the class.
-
- c: instance of class Control
- Contains all the parameters of control files, which are e.g.:
- DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
- STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
- LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
- OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
- ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
- MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
- DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
-
- For more information about format and content of the parameter
- see documentation.
-
- fluxes: boolean, optional
- Decides if a the flux parameter settings are stored or
- the rest of the parameter list.
- Default value is False.
-
- @Return:
- <nothing>
- '''
-
- # different mars types for retrieving reanalysis data for flexpart
- self.types = dict()
- try:
- if c.maxstep > len(c.type): # Pure forecast mode
- c.type = [c.type[1]]
- c.step = ['{:0>3}'.format(int(c.step[0]))]
- c.time = [c.time[0]]
- for i in range(1, c.maxstep + 1):
- c.type.append(c.type[0])
- c.step.append('{:0>3}'.format(i))
- c.time.append(c.time[0])
- except:
- pass
-
- self.inputdir = c.inputdir
- self.basetime = c.basetime
- self.dtime = c.dtime
- self.mars = {}
- i = 0
- j = 0
- if fluxes is True and c.maxstep < 24:
- # only FC with start times at 00/12 (without 00/12)
- self.types[c.type[1]] = {'times': '00/12',
- 'steps': '{}/to/12/by/{}'.format(c.dtime,
- c.dtime)}
- i = 1
- for k in [0, 12]:
- for j in range(int(c.dtime), 13, int(c.dtime)):
- self.mars['{:0>3}'.format(i * int(c.dtime))] = \
- [c.type[1], '{:0>2}'.format(k), '{:0>3}'.format(j)]
- i += 1
- else:
- for ty, st, ti in zip(c.type, c.step, c.time):
- btlist = range(24)
- if c.basetime == '12':
- btlist = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
- if c.basetime == '00':
- btlist = [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 0]
-
- if mod(i, int(c.dtime)) == 0 and \
- (i in btlist or c.maxstep > 24):
-
- if ty not in self.types.keys():
- self.types[ty] = {'times': '', 'steps': ''}
-
- if ti not in self.types[ty]['times']:
- if len(self.types[ty]['times']) > 0:
- self.types[ty]['times'] += '/'
- self.types[ty]['times'] += ti
-
- if st not in self.types[ty]['steps']:
- if len(self.types[ty]['steps']) > 0:
- self.types[ty]['steps'] += '/'
- self.types[ty]['steps'] += st
-
- self.mars['{:0>3}'.format(j)] = [ty,
- '{:0>2}'.format(int(ti)),
- '{:0>3}'.format(int(st))]
- j += int(c.dtime)
-
- i += 1
-
- # Different grids need different retrievals
- # SH = Spherical Harmonics, GG = Gaussian Grid,
- # OG = Output Grid, ML = MultiLevel, SL = SingleLevel
- self.params = {'SH__ML': '', 'SH__SL': '',
- 'GG__ML': '', 'GG__SL': '',
- 'OG__ML': '', 'OG__SL': '',
- 'OG_OROLSM_SL': '', 'OG_acc_SL': ''}
- self.marsclass = c.marsclass
- self.stream = c.stream
- self.number = c.number
- self.resol = c.resol
- if 'N' in c.grid: # Gaussian output grid
- self.grid = c.grid
- self.area = 'G'
- else:
- self.grid = '{}/{}'.format(int(c.grid) / 1000., int(c.grid) / 1000.)
- self.area = '{}/{}/{}/{}'.format(int(c.upper) / 1000.,
- int(c.left) / 1000.,
- int(c.lower) / 1000.,
- int(c.right) / 1000.)
-
- self.accuracy = c.accuracy
- self.level = c.level
- try:
- self.levelist = c.levelist
- except:
- self.levelist = '1/to/' + c.level
-
- self.glevelist = '1/to/' + c.level
-
- try:
- self.gaussian = c.gaussian
- except:
- self.gaussian = ''
-
- try:
- self.dataset = c.dataset
- except:
- self.dataset = ''
-
- try:
- self.expver = c.expver
- except:
- self.expver = '1'
-
- try:
- self.number = c.number
- except:
- self.number = '0'
-
- self.outputfilelist = []
-
-
- # Now comes the nasty part that deals with the different
- # scenarios we have:
- # 1) Calculation of etadot on
- # a) Gaussian grid
- # b) Output grid
- # c) Output grid using parameter 77 retrieved from MARS
- # 3) Calculation/Retrieval of omega
- # 4) Download also data for WRF
-
- if fluxes is False:
- self.params['SH__SL'] = ['LNSP', 'ML', '1', 'OFF']
- # self.params['OG__SL'] = ["SD/MSL/TCC/10U/10V/2T/2D/129/172",
- # 'SFC','1',self.grid]
- self.params['OG__SL'] = ["141/151/164/165/166/167/168/129/172", \
- 'SFC', '1', self.grid]
- self.params['OG_OROLSM__SL'] = ["160/27/28/173", \
- 'SFC', '1', self.grid]
-
- if len(c.addpar) > 0:
- if c.addpar[0] == '/':
- c.addpar = c.addpar[1:]
- self.params['OG__SL'][0] += '/' + '/'.join(c.addpar)
- self.params['OG__ML'] = ['T/Q', 'ML', self.levelist, self.grid]
-
- if c.gauss == '0' and c.eta == '1':
- # the simplest case
- self.params['OG__ML'][0] += '/U/V/77'
- elif c.gauss == '0' and c.eta == '0':
- # this is not recommended (inaccurate)
- self.params['OG__ML'][0] += '/U/V'
- elif c.gauss == '1' and c.eta == '0':
- # this is needed for data before 2008, or for reanalysis data
- self.params['GG__SL'] = ['Q', 'ML', '1', \
- '{}'.format((int(self.resol) + 1) / 2)]
- self.params['SH__ML'] = ['U/V/D', 'ML', self.glevelist, 'OFF']
- else:
- print('Warning: This is a very costly parameter combination, \
- use only for debugging!')
- self.params['GG__SL'] = ['Q', 'ML', '1', \
- '{}'.format((int(self.resol) + 1) / 2)]
- self.params['GG__ML'] = ['U/V/D/77', 'ML', self.glevelist, \
- '{}'.format((int(self.resol) + 1) / 2)]
-
- if c.omega == '1':
- self.params['OG__ML'][0] += '/W'
-
- try:
- # add cloud water content if necessary
- if c.cwc == '1':
- self.params['OG__ML'][0] += '/CLWC/CIWC'
- except:
- pass
-
- try:
- # add vorticity and geopotential height for WRF if necessary
- if c.wrf == '1':
- self.params['OG__ML'][0] += '/Z/VO'
- if '/D' not in self.params['OG__ML'][0]:
- self.params['OG__ML'][0] += '/D'
- #wrf_sfc = 'sp/msl/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/stl1/ /
- # stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4'.upper()
- wrf_sfc = '134/235/167/165/166/168/129/172/34/31/141/ \
- 139/170/183/236/39/40/41/42'.upper()
- lwrt_sfc = wrf_sfc.split('/')
- for par in lwrt_sfc:
- if par not in self.params['OG__SL'][0]:
- self.params['OG__SL'][0] += '/' + par
- except:
- pass
- else:
- self.params['OG_acc_SL'] = ["LSP/CP/SSHF/EWSS/NSSS/SSR", \
- 'SFC', '1', self.grid]
-
- # if needed, add additional WRF specific parameters here
-
- return
-
-
- def write_namelist(self, c, filename): #done
- '''
- @Description:
- Creates a namelist file in the temporary directory and writes
- the following values to it: maxl, maxb, mlevel,
- mlevelist, mnauf, metapar, rlo0, rlo1, rla0, rla1,
- momega, momegadiff, mgauss, msmooth, meta, metadiff, mdpdeta
-
- @Input:
- self: instance of ECFlexpart
- The current object of the class.
-
- c: instance of class Control
- Contains all the parameters of control files, which are e.g.:
- DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
- STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
- LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
- OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
- ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
- MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
- DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
-
- For more information about format and content of the parameter
- see documentation.
-
- filename: string
- Name of the namelist file.
-
- @Return:
- <nothing>
- '''
-
- self.inputdir = c.inputdir
- area = asarray(self.area.split('/')).astype(float)
- grid = asarray(self.grid.split('/')).astype(float)
-
- if area[1] > area[3]:
- area[1] -= 360
- zyk = abs((area[3] - area[1] - 360.) + grid[1]) < 1.e-6
- maxl = int((area[3] - area[1]) / grid[1]) + 1
- maxb = int((area[0] - area[2]) / grid[0]) + 1
-
- with open(self.inputdir + '/' + filename, 'w') as f:
- f.write('&NAMGEN\n')
- f.write(',\n '.join(['maxl = ' + str(maxl), 'maxb = ' + str(maxb),
- 'mlevel = ' + self.level,
- 'mlevelist = ' + '"' + self.levelist + '"',
- 'mnauf = ' + self.resol, 'metapar = ' + '77',
- 'rlo0 = ' + str(area[1]), 'rlo1 = ' + str(area[3]),
- 'rla0 = ' + str(area[2]), 'rla1 = ' + str(area[0]),
- 'momega = ' + c.omega, 'momegadiff = ' + c.omegadiff,
- 'mgauss = ' + c.gauss, 'msmooth = ' + c.smooth,
- 'meta = ' + c.eta, 'metadiff = ' + c.etadiff,
- 'mdpdeta = ' + c.dpdeta]))
-
- f.write('\n/\n')
-
- return
-
- def retrieve(self, server, dates, times, inputdir=''):
- '''
- @Description:
-
-
- @Input:
- self: instance of ECFlexpart
-
- server: instance of ECMWFService
-
- dates:
-
- times:
-
- inputdir: string, optional
- Default string is empty ('').
-
- @Return:
- <nothing>
- '''
- self.dates = dates
- self.server = server
-
- if inputdir == "":
- self.inputdir = '.'
- else:
- self.inputdir = inputdir
-
- # Retrieve Q not for using Q but as a template for a reduced gaussian
- # grid one date and time is enough
- # Take analysis at 00
- qdate = self.dates
- idx = qdate.find("/")
- if (idx > 0):
- qdate = self.dates[:idx]
-
- #QG = MARSretrieval(self.server, dataset = self.dataset, marsclass = self.marsclass, stream = self.stream, type = "an", levtype = "ML", levelist = "1",
- #gaussian = "reduced",grid = '{}'.format((int(self.resol)+1)/2), resol = self.resol,accuracy = self.accuracy,target = self.inputdir+"/"+"QG.grb",
- #date = qdate, time = "00",expver = self.expver, param = "133.128")
- #QG.displayInfo()
- #QG.dataRetrieve()
-
- oro = False
- for ftype in self.types:
- for pk, pv in self.params.iteritems():
- if isinstance(pv, str): # or pk != 'GG__SL' :
- continue
- mftype = ''+ftype
- mftime = self.types[ftype]['times']
- mfstep = self.types[ftype]['steps']
- mfdate = self.dates
- mfstream = self.stream
- mftarget = self.inputdir + "/" + ftype + pk + '.' + \
- self.dates.split('/')[0] + '.' + str(os.getppid()) +\
- '.' + str(os.getpid()) + ".grb"
- if pk == 'OG__SL':
- pass
- if pk == 'OG_OROLSM__SL':
- if oro is False:
- mfstream = 'OPER'
- mftype = 'AN'
- mftime = '00'
- mfstep = '000'
- mfdate = self.dates.split('/')[0]
- mftarget = self.inputdir + "/" + pk + '.' + mfdate + \
- '.' + str(os.getppid()) + '.' + \
- str(os.getpid()) + ".grb"
- oro = True
- else:
- continue
-
- if pk == 'GG__SL' and pv[0] == 'Q':
- area = ""
- gaussian = 'reduced'
- else:
- area = self.area
- gaussian = self.gaussian
-
- if self.basetime is None:
- MR = MARSretrieval(self.server,
- marsclass=self.marsclass, stream=mfstream,
- type=mftype, levtype=pv[1], levelist=pv[2],
- resol=self.resol, gaussian=gaussian,
- accuracy=self.accuracy, grid=pv[3],
- target=mftarget, area=area, date=mfdate,
- time=mftime, number=self.number, step=mfstep,
- expver=self.expver, param=pv[0])
-
- MR.displayInfo()
- MR.dataRetrieve()
- # The whole else section is only necessary for operational scripts.
- # It could be removed
- else:
- # check if mars job requests fields beyond basetime.
- # If yes eliminate those fields since they may not
- # be accessible with user's credentials
- sm1 = -1
- if 'by' in mfstep:
- sm1 = 2
- tm1 = -1
- if 'by' in mftime:
- tm1 = 2
- maxtime = datetime.datetime.strptime(
- mfdate.split('/')[-1] + mftime.split('/')[tm1],
- '%Y%m%d%H') + datetime.timedelta(
- hours=int(mfstep.split('/')[sm1]))
-
- elimit = datetime.datetime.strptime(
- mfdate.split('/')[-1] +
- self.basetime, '%Y%m%d%H')
-
- if self.basetime == '12':
- if 'acc' in pk:
-
- # Strategy: if maxtime-elimit> = 24h reduce date by 1,
- # if 12h< = maxtime-elimit<12h reduce time for last date
- # if maxtime-elimit<12h reduce step for last time
- # A split of the MARS job into 2 is likely necessary.
- maxtime = elimit-datetime.timedelta(hours=24)
- mfdate = '/'.join(('/'.join(mfdate.split('/')[:-1]),
- datetime.datetime.strftime(
- maxtime, '%Y%m%d')))
-
- MR = MARSretrieval(self.server,
- marsclass=self.marsclass,
- stream=self.stream, type=mftype,
- levtype=pv[1], levelist=pv[2],
- resol=self.resol, gaussian=gaussian,
- accuracy=self.accuracy, grid=pv[3],
- target=mftarget, area=area,
- date=mfdate, time=mftime,
- number=self.number, step=mfstep,
- expver=self.expver, param=pv[0])
-
- MR.displayInfo()
- MR.dataRetrieve()
-
- maxtime = elimit - datetime.timedelta(hours=12)
- mfdate = datetime.datetime.strftime(maxtime,
- '%Y%m%d')
- mftime = '00'
- mftarget = self.inputdir + "/" + ftype + pk + \
- '.' + mfdate + '.' + str(os.getppid()) +\
- '.' + str(os.getpid()) + ".grb"
-
- MR = MARSretrieval(self.server,
- marsclass=self.marsclass,
- stream=self.stream, type=mftype,
- levtype=pv[1], levelist=pv[2],
- resol=self.resol, gaussian=gaussian,
- accuracy=self.accuracy, grid=pv[3],
- target=mftarget, area=area,
- date=mfdate, time=mftime,
- number=self.number, step=mfstep,
- expver=self.expver, param=pv[0])
-
- MR.displayInfo()
- MR.dataRetrieve()
- else:
- MR = MARSretrieval(self.server,
- marsclass=self.marsclass,
- stream=self.stream, type=mftype,
- levtype=pv[1], levelist=pv[2],
- resol=self.resol, gaussian=gaussian,
- accuracy=self.accuracy, grid=pv[3],
- target=mftarget, area=area,
- date=mfdate, time=mftime,
- number=self.number, step=mfstep,
- expver=self.expver, param=pv[0])
-
- MR.displayInfo()
- MR.dataRetrieve()
- else:
- maxtime = elimit - datetime.timedelta(hours=24)
- mfdate = datetime.datetime.strftime(maxtime,'%Y%m%d')
-
- mftimesave = ''.join(mftime)
-
- if '/' in mftime:
- times = mftime.split('/')
- while ((int(times[0]) +
- int(mfstep.split('/')[0]) <= 12) and
- (pk != 'OG_OROLSM__SL') and 'acc' not in pk):
- times = times[1:]
- if len(times) > 1:
- mftime = '/'.join(times)
- else:
- mftime = times[0]
-
- MR = MARSretrieval(self.server,
- marsclass=self.marsclass,
- stream=self.stream, type=mftype,
- levtype=pv[1], levelist=pv[2],
- resol=self.resol, gaussian=gaussian,
- accuracy=self.accuracy, grid=pv[3],
- target=mftarget, area=area,
- date=mfdate, time=mftime,
- number=self.number, step=mfstep,
- expver=self.expver, param=pv[0])
-
- MR.displayInfo()
- MR.dataRetrieve()
-
- if (int(mftimesave.split('/')[0]) == 0 and
- int(mfstep.split('/')[0]) == 0 and
- pk != 'OG_OROLSM__SL'):
- mfdate = datetime.datetime.strftime(elimit,'%Y%m%d')
- mftime = '00'
- mfstep = '000'
- mftarget = self.inputdir + "/" + ftype + pk + \
- '.' + mfdate + '.' + str(os.getppid()) +\
- '.' + str(os.getpid()) + ".grb"
-
- MR = MARSretrieval(self.server,
- marsclass=self.marsclass,
- stream=self.stream, type=mftype,
- levtype=pv[1], levelist=pv[2],
- resol=self.resol, gaussian=gaussian,
- accuracy=self.accuracy, grid=pv[3],
- target=mftarget, area=area,
- date=mfdate, time=mftime,
- number=self.number, step=mfstep,
- expver=self.expver, param=pv[0])
-
- MR.displayInfo()
- MR.dataRetrieve()
-
- print("MARS retrieve done... ")
-
- return
-
-
- def process_output(self, c): #done
- '''
- @Description:
- The grib files are postprocessed depending on selection in
- control file. The following modifications might be done if
- properly switched in control file:
- GRIB2 - Conversion to GRIB2
- ECTRANS - Transfer of files to gateway server
- ECSTORAGE - Storage at ECMWF server
- GRIB2FLEXPART - Conversion of GRIB files to FLEXPART binary format
-
- @Input:
- self: instance of ECFlexpart
- The current object of the class.
-
- c: instance of class Control
- Contains all the parameters of control files, which are e.g.:
- DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
- STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
- LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
- OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
- ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
- MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
- DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
-
- For more information about format and content of the parameter
- see documentation.
-
- @Return:
- <nothing>
-
- '''
-
- print('Postprocessing:\n Format: {}\n'.format(c.format))
-
- if c.ecapi is False:
- print('ecstorage: {}\n ecfsdir: {}\n'.
- format(c.ecstorage, c.ecfsdir))
- if not hasattr(c, 'gateway'):
- c.gateway = os.getenv('GATEWAY')
- if not hasattr(c, 'destination'):
- c.destination = os.getenv('DESTINATION')
- print('ectrans: {}\n gateway: {}\n destination: {}\n '
- .format(c.ectrans, c.gateway, c.destination))
-
- print('Output filelist: \n', self.outputfilelist)
-
- if c.format.lower() == 'grib2':
- for ofile in self.outputfilelist:
- p = subprocess.check_call(['grib_set', '-s', 'edition=2, \
- productDefinitionTemplateNumber=8',
- ofile, ofile + '_2'])
- p = subprocess.check_call(['mv', ofile + '_2', ofile])
-
- if int(c.ectrans) == 1 and c.ecapi is False:
- for ofile in self.outputfilelist:
- p = subprocess.check_call(['ectrans', '-overwrite', '-gateway',
- c.gateway, '-remote', c.destination,
- '-source', ofile])
- print('ectrans:', p)
-
- if int(c.ecstorage) == 1 and c.ecapi is False:
- for ofile in self.outputfilelist:
- p = subprocess.check_call(['ecp', '-o', ofile,
- os.path.expandvars(c.ecfsdir)])
-
- if c.outputdir != c.inputdir:
- for ofile in self.outputfilelist:
- p = subprocess.check_call(['mv', ofile, c.outputdir])
-
- # prepare environment for the grib2flexpart run
- # to convert grib to flexpart binary
- if c.grib2flexpart == '1':
-
- # generate AVAILABLE file
- # Example of AVAILABLE file data
- # 20131107 000000 EN13110700 ON DISC
- clist = []
- for ofile in self.outputfilelist:
- fname = ofile.split('/')
- if '.' in fname[-1]:
- l = fname[-1].split('.')
- timestamp = datetime.datetime.strptime(l[0][-6:] + l[1],
- '%y%m%d%H')
- timestamp += datetime.timedelta(hours=int(l[2]))
- cdate = datetime.datetime.strftime(timestamp, '%Y%m%d')
- chms = datetime.datetime.strftime(timestamp, '%H%M%S')
- else:
- cdate = '20' + fname[-1][-8:-2]
- chms = fname[-1][-2:] + '0000'
- clist.append(cdate + ' ' + chms + ' '*6 +
- fname[-1] + ' '*14 + 'ON DISC')
- clist.sort()
- with open(c.outputdir + '/' + 'AVAILABLE', 'w') as f:
- f.write('\n'.join(clist) + '\n')
-
- # generate pathnames file
- pwd = os.path.abspath(c.outputdir)
- with open(pwd + '/pathnames','w') as f:
- f.write(pwd + '/Options/\n')
- f.write(pwd + '/\n')
- f.write(pwd + '/\n')
- f.write(pwd + '/AVAILABLE\n')
- f.write(' = == = == = == = == = == == = \n')
-
- # create Options dir if necessary
- if not os.path.exists(pwd + '/Options'):
- os.makedirs(pwd+'/Options')
-
- # read template COMMAND file
- with open(os.path.expandvars(
- os.path.expanduser(c.flexpart_root_scripts)) +
- '/../Options/COMMAND', 'r') as f:
- lflist = f.read().split('\n')
-
- # find index of list where to put in the
- # date and time information
- # usually after the LDIRECT parameter
- i = 0
- for l in lflist:
- if 'LDIRECT' in l.upper():
- break
- i += 1
-
-# clist.sort()
- # insert the date and time information of run star and end
- # into the list of lines of COMMAND file
- lflist = lflist[:i+1] + \
- [clist[0][:16], clist[-1][:16]] + \
- lflist[i+3:]
-
- # write the new COMMAND file
- with open(pwd + '/Options/COMMAND', 'w') as g:
- g.write('\n'.join(lflist) + '\n')
-
- # change to outputdir and start the
- # grib2flexpart run
- # afterwards switch back to the working dir
- os.chdir(c.outputdir)
- p = subprocess.check_call([os.path.expandvars(
- os.path.expanduser(c.flexpart_root_scripts)) +
- '/../FLEXPART_PROGRAM/grib2flexpart',
- 'useAvailable', '.'])
- os.chdir(pwd)
-
- return
-
- def create(self, inputfiles, c): #done
- '''
- @Description:
- This method is based on the ECMWF example index.py
- https://software.ecmwf.int/wiki/display/GRIB/index.py
-
- An index file will be created which depends on the combination
- of "date", "time" and "stepRange" values. This is used to iterate
- over all messages in the grib files passed through the parameter
- "inputfiles" to seperate specific parameters into fort.* files.
- Afterwards the FORTRAN program Convert2 is called to convert
- the data fields all to the same grid and put them in one file
- per day.
-
- @Input:
- self: instance of ECFlexpart
- The current object of the class.
-
- inputfiles: instance of UIOFiles
- Contains a list of files.
-
- c: instance of class Control
- Contains all the parameters of control files, which are e.g.:
- DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
- STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
- LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
- OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
- ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
- MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
- DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
-
- For more information about format and content of the parameter
- see documentation.
-
- @Return:
- <nothing>
- '''
-
- table128 = init128(c.ecmwfdatadir +
- '/grib_templates/ecmwf_grib1_table_128')
- wrfpars = toparamId('sp/mslp/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/\
- stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4',
- table128)
-
- index_keys = ["date", "time", "step"]
- indexfile = c.inputdir + "/date_time_stepRange.idx"
- silentremove(indexfile)
- grib = GribTools(inputfiles.files)
- # creates new index file
- iid = grib.index(index_keys=index_keys, index_file=indexfile)
-
- # read values of index keys
- index_vals = []
- for key in index_keys:
- index_vals.append(grib_index_get(iid, key))
- print(index_vals[-1])
- # index_vals looks for example like:
- # index_vals[0]: ('20171106', '20171107', '20171108') ; date
- # index_vals[1]: ('0', '1200', '1800', '600') ; time
- # index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange
-
-
- for prod in product(*index_vals):
- # e.g. prod = ('20170505', '0', '12')
- # ( date ,time, step)
- # per date e.g. time = 0, 600, 1200, 1800
- # per time e.g. step = 0, 3, 6, 9, 12
- for i in range(len(index_keys)):
- grib_index_select(iid, index_keys[i], prod[i])
-
- gid = grib_new_from_index(iid)
- # do convert2 program if gid at this time is not None,
- # therefore save in hid
- hid = gid
- if gid is not None:
- cdate = str(grib_get(gid, 'date'))
- time = grib_get(gid, 'time')
- type = grib_get(gid, 'type')
- step = grib_get(gid, 'step')
- # create correct timestamp from the three time informations
- # date, time, step
- timestamp = datetime.datetime.strptime(
- cdate + '{:0>2}'.format(time/100), '%Y%m%d%H')
- timestamp += datetime.timedelta(hours=int(step))
-
- cdateH = datetime.datetime.strftime(timestamp, '%Y%m%d%H')
- chms = datetime.datetime.strftime(timestamp, '%H%M%S')
-
- if c.basetime is not None:
- slimit = datetime.datetime.strptime(
- c.start_date + '00', '%Y%m%d%H')
- bt = '23'
- if c.basetime == '00':
- bt = '00'
- slimit = datetime.datetime.strptime(
- c.end_date + bt, '%Y%m%d%H') - \
- datetime.timedelta(hours=12-int(c.dtime))
- if c.basetime == '12':
- bt = '12'
- slimit = datetime.datetime.strptime(
- c.end_date + bt, '%Y%m%d%H') - \
- datetime.timedelta(hours=12-int(c.dtime))
-
- elimit = datetime.datetime.strptime(
- c.end_date + bt, '%Y%m%d%H')
-
- if timestamp < slimit or timestamp > elimit:
- continue
-
- try:
- if c.wrf == '1':
- if 'olddate' not in locals():
- fwrf = open(c.outputdir + '/WRF' + cdate +
- '.{:0>2}'.format(time) + '.000.grb2', 'w')
- olddate = cdate[:]
- else:
- if cdate != olddate:
- fwrf = open(c.outputdir + '/WRF' + cdate +
- '.{:0>2}'.format(time) + '.000.grb2',
- 'w')
- olddate = cdate[:]
- except AttributeError:
- pass
-
- # delete old fort.* files and open them newly
- fdict = {'10':None, '11':None, '12':None, '13':None, '16':None,
- '17':None, '19':None, '21':None, '22':None, '20':None}
- #for f in fdict.keys():
- # silentremove(c.inputdir + "/fort." + f)
- for k, f in fdict.iteritems():
- silentremove(c.inputdir + "/fort." + k)
- fdict[k] = open(c.inputdir + '/fort.' + k, 'w')
-
- savedfields = []
- while 1:
- if gid is None:
- break
- paramId = grib_get(gid, 'paramId')
- gridtype = grib_get(gid, 'gridType')
- datatype = grib_get(gid, 'dataType')
- levtype = grib_get(gid, 'typeOfLevel')
- if paramId == 133 and gridtype == 'reduced_gg':
- # Relative humidity (Q.grb) is used as a template only
- # so we need the first we "meet"
- with open(c.inputdir + '/fort.18', 'w') as fout:
- grib_write(gid, fout)
- elif paramId == 131 or paramId == 132:
- grib_write(gid, fdict['10'])
- elif paramId == 130:
- grib_write(gid, fdict['11'])
- elif paramId == 133 and gridtype != 'reduced_gg':
- grib_write(gid, fdict['17'])
- elif paramId == 152:
- grib_write(gid, fdict['12'])
- elif paramId == 155 and gridtype == 'sh':
- grib_write(gid, fdict['13'])
- elif paramId in [129, 138, 155] and levtype == 'hybrid' \
- and c.wrf == '1':
- pass
- elif paramId == 246 or paramId == 247:
- # cloud liquid water and ice
- if paramId == 246:
- clwc = grib_get_values(gid)
- else:
- clwc += grib_get_values(gid)
- grib_set_values(gid, clwc)
- grib_set(gid, 'paramId', 201031)
- grib_write(gid, fdict['22'])
- elif paramId == 135:
- grib_write(gid, fdict['19'])
- elif paramId == 77:
- grib_write(gid, fdict['21'])
- else:
- if paramId not in savedfields:
- grib_write(gid, fdict['16'])
- savedfields.append(paramId)
- else:
- print('duplicate ' + str(paramId) + ' not written')
-
- try:
- if c.wrf == '1':
-# die if abfrage scheint ueberfluessig da eh das gleihce ausgefuehrt wird
- if levtype == 'hybrid':
- if paramId in [129, 130, 131, 132, 133, 138, 155]:
- grib_write(gid, fwrf)
- else:
- if paramId in wrfpars:
- grib_write(gid, fwrf)
- except AttributeError:
- pass
-
- grib_release(gid)
- gid = grib_new_from_index(iid)
-
- for f in fdict.values():
- f.close()
-
- # call for CONVERT2
-# AUSLAGERN IN EIGENE FUNKTION
-
- if hid is not None:
- pwd = os.getcwd()
- os.chdir(c.inputdir)
- if os.stat('fort.21').st_size == 0 and int(c.eta) == 1:
- print('Parameter 77 (etadot) is missing, most likely it is \
- not available for this type or date/time\n')
- print('Check parameters CLASS, TYPE, STREAM, START_DATE\n')
- myerror(c, 'fort.21 is empty while parameter eta is set \
- to 1 in CONTROL file')
-
- p = subprocess.check_call([os.path.expandvars(
- os.path.expanduser(c.exedir)) + '/CONVERT2'], shell=True)
- os.chdir(pwd)
- # create the corresponding output file fort.15
- # (generated by CONVERT2)
- # + fort.16 (paramId 167 and paramId 168)
- fnout = c.inputdir + '/' + c.prefix
- if c.maxstep > 12:
- suffix = cdate[2:8] + '.{:0>2}'.format(time/100) + \
- '.{:0>3}'.format(step)
- else:
- suffix = cdateH[2:10]
-
- fnout += suffix
- print("outputfile = " + fnout)
- self.outputfilelist.append(fnout) # needed for final processing
- fout = open(fnout, 'wb')
- shutil.copyfileobj(open(c.inputdir + '/fort.15', 'rb'), fout)
- if c.cwc == '1':
- shutil.copyfileobj(open(c.inputdir + '/fort.22', 'rb'), fout)
- shutil.copyfileobj(open(c.inputdir + '/flux' + cdate[0:2] +
- suffix, 'rb'), fout)
- shutil.copyfileobj(open(c.inputdir + '/fort.16', 'rb'), fout)
- orolsm = glob.glob(c.inputdir +
- '/OG_OROLSM__SL.*.' + c.ppid + '*')[0]
- shutil.copyfileobj(open(orolsm, 'rb'), fout)
- fout.close()
- if c.omega == '1':
- fnout = c.outputdir + '/OMEGA'
- fout = open(fnout, 'wb')
- shutil.copyfileobj(open(c.inputdir + '/fort.25', 'rb'), fout)
-
- try:
- if c.wrf == '1':
- fwrf.close()
- except:
- pass
-
- grib_index_release(iid)
-
- return
-
-
- def deacc_fluxes(self, inputfiles, c):
- '''
- @Description:
-
-
- @Input:
- self: instance of ECFlexpart
- The current object of the class.
-
- inputfiles: instance of UIOFiles
- Contains a list of files.
-
- c: instance of class Control
- Contains all the parameters of control files, which are e.g.:
- DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
- STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
- LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
- OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
- ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
- MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
- DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
-
- For more information about format and content of the parameter
- see documentation.
-
- @Return:
- <nothing>
- '''
-
- table128 = init128(c.ecmwfdatadir +
- '/grib_templates/ecmwf_grib1_table_128')
- pars = toparamId(self.params['OG_acc_SL'][0], table128)
-
- index_keys = ["date", "time", "step"]
- indexfile = c.inputdir + "/date_time_stepRange.idx"
- silentremove(indexfile)
- grib = GribTools(inputfiles.files)
- # creates new index file
- iid = grib.index(index_keys=index_keys, index_file=indexfile)
-
- # read values of index keys
- index_vals = []
- for key in index_keys:
- key_vals = grib_index_get(iid,key)
- print(key_vals)
- # have to sort the steps for disaggregation,
- # therefore convert to int first
- if key == 'step':
- key_vals = [int(k) for k in key_vals]
- key_vals.sort()
- key_vals = [str(k) for k in key_vals]
- index_vals.append(key_vals)
- # index_vals looks for example like:
- # index_vals[0]: ('20171106', '20171107', '20171108') ; date
- # index_vals[1]: ('0', '1200', '1800', '600') ; time
- # index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange
-
- valsdict = {}
- svalsdict = {}
- stepsdict = {}
- for p in pars:
- valsdict[str(p)] = []
- svalsdict[str(p)] = []
- stepsdict[str(p)] = []
- # ab hier eien Einrückung zurück!!!!
- for prod in product(*index_vals):
- # e.g. prod = ('20170505', '0', '12')
- # ( date ,time, step)
- # per date e.g. time = 0, 600, 1200, 1800
- # per time e.g. step = 0, 3, 6, 9, 12
- for i in range(len(index_keys)):
- grib_index_select(iid, index_keys[i], prod[i])
-
- gid = grib_new_from_index(iid)
- # do convert2 program if gid at this time is not None,
- # therefore save in hid
- hid = gid
- if gid is not None:
- cdate = str(grib_get(gid, 'date'))
- time = grib_get(gid, 'time')
- type = grib_get(gid, 'type')
- step = grib_get(gid, 'step')
- # date+time+step-2*dtime
- #(since interpolated value valid for step-2*dtime)
- sdate = datetime.datetime(year=cdate / 10000,
- month=mod(cdate, 10000) / 100,
- day=mod(cdate, 100),
- hour=time / 100)
- fdate = sdate + datetime.timedelta(
- hours=step - 2 * int(c.dtime))
- sdates = sdate + datetime.timedelta(hours=step)
- else:
- break
-
- if c.maxstep > 12:
- fnout = c.inputdir + '/flux' + \
- sdate.strftime('%Y%m%d') + '.{:0>2}'.format(time/100) + \
- '.{:0>3}'.format(step-2*int(c.dtime))
- gnout = c.inputdir + '/flux' + \
- sdate.strftime('%Y%m%d')+'.{:0>2}'.format(time/100) + \
- '.{:0>3}'.format(step-int(c.dtime))
- hnout = c.inputdir + '/flux' + \
- sdate.strftime('%Y%m%d')+'.{:0>2}'.format(time/100) + \
- '.{:0>3}'.format(step)
- g = open(gnout, 'w')
- h = open(hnout, 'w')
- else:
- fnout = c.inputdir + '/flux' + fdate.strftime('%Y%m%d%H')
- gnout = c.inputdir + '/flux' + (fdate+datetime.timedelta(
- hours = int(c.dtime))).strftime('%Y%m%d%H')
- hnout = c.inputdir + '/flux' + sdates.strftime('%Y%m%d%H')
- g = open(gnout, 'w')
- h = open(hnout, 'w')
- print("outputfile = " + fnout)
- f = open(fnout, 'w')
-
- while 1:
- if gid is None:
- break
- cparamId = str(grib_get(gid, 'paramId'))
- step = grib_get(gid, 'step')
- atime = grib_get(gid, 'time')
- ni = grib_get(gid, 'Ni')
- nj = grib_get(gid, 'Nj')
- if cparamId in valsdict.keys():
- values = grib_get_values(gid)
- vdp = valsdict[cparamId]
- svdp = svalsdict[cparamId]
- sd = stepsdict[cparamId]
-
- if cparamId == '142' or cparamId == '143':
- fak = 1./1000.
- else:
- fak = 3600.
-
- values = (reshape(values, (nj, ni))).flatten() / fak
- vdp.append(values[:]) # save the accumulated values
- if step <= int(c.dtime):
- svdp.append(values[:] / int(c.dtime))
- else: # deaccumulate values
- svdp.append((vdp[-1] - vdp[-2]) / int(c.dtime))
-
- print(cparamId, atime, step, len(values),
- values[0], std(values))
- # save the 1/3-hourly or specific values
- # svdp.append(values[:])
- sd.append(step)
- # len(svdp) correspond to the time
- if len(svdp) >= 3:
- if len(svdp) > 3:
- if cparamId == '142' or cparamId == '143':
- values = darain(svdp)
- else:
- values = dapoly(svdp)
-
- if not (step == c.maxstep and c.maxstep > 12 \
- or sdates == elimit):
- vdp.pop(0)
- svdp.pop(0)
- else:
- if c.maxstep > 12:
- values = svdp[1]
- else:
- values = svdp[0]
-
- grib_set_values(gid, values)
- if c.maxstep > 12:
- grib_set(gid, 'step', max(0, step-2*int(c.dtime)))
- else:
- grib_set(gid, 'step', 0)
- grib_set(gid, 'time', fdate.hour*100)
- grib_set(gid, 'date', fdate.year*10000 +
- fdate.month*100+fdate.day)
- grib_write(gid, f)
-
- if c.basetime is not None:
- elimit = datetime.datetime.strptime(c.end_date +
- c.basetime,
- '%Y%m%d%H')
- else:
- elimit = sdate + datetime.timedelta(2*int(c.dtime))
-
- # squeeze out information of last two steps contained
- # in svdp
- # if step+int(c.dtime) == c.maxstep and c.maxstep>12
- # or sdates+datetime.timedelta(hours = int(c.dtime))
- # >= elimit:
- # Note that svdp[0] has not been popped in this case
-
- if (step == c.maxstep and c.maxstep > 12
- or sdates == elimit):
- values = svdp[3]
- grib_set_values(gid, values)
- grib_set(gid, 'step', 0)
- truedatetime = fdate + datetime.timedelta(
- hours=2*int(c.dtime))
- grib_set(gid, 'time', truedatetime.hour * 100)
- grib_set(gid, 'date', truedatetime.year * 10000 +
- truedatetime.month * 100 +
- truedatetime.day)
- grib_write(gid, h)
-
- #values = (svdp[1]+svdp[2])/2.
- if cparamId == '142' or cparamId == '143':
- values = darain(list(reversed(svdp)))
- else:
- values = dapoly(list(reversed(svdp)))
-
- grib_set(gid, 'step',0)
- truedatetime = fdate + datetime.timedelta(
- hours=int(c.dtime))
- grib_set(gid, 'time', truedatetime.hour * 100)
- grib_set(gid, 'date', truedatetime.year * 10000 +
- truedatetime.month * 100 +
- truedatetime.day)
- grib_set_values(gid, values)
- grib_write(gid, g)
-
- grib_release(gid)
-
- gid = grib_new_from_index(iid)
-
- f.close()
- g.close()
- h.close()
-
- grib_index_release(iid)
-
- return
-
diff --git a/python/GribTools.py b/python/GribTools.py
index 29b7626..3029352 100644
--- a/python/GribTools.py
+++ b/python/GribTools.py
@@ -5,7 +5,8 @@
#AP
# - GribTools name möglicherweise etwas verwirrend.
# - change self.filename in self.filenames!!!
-# -
+# - add file description
+# - bis auf --init-- und index wird keine Funktion verwendet!?
#************************************************************************
"""
@Author: Anne Fouilloux (University of Oslo)
@@ -19,7 +20,7 @@
- changed some naming
@License:
- (C) Copyright 2014 UIO.
+ (C) Copyright 2014-2018.
This software is licensed under the terms of the Apache Licence Version 2.0
which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
@@ -27,7 +28,6 @@
@Requirements:
- A standard python 2.6 or 2.7 installation
- dateutils
- - matplotlib (optional, for debugging)
- ECMWF specific packages, all available from https://software.ecmwf.int/
ECMWF WebMARS, gribAPI with python enabled, emoslib and
ecaccess web toolkit
@@ -35,17 +35,15 @@
@Description:
Further documentation may be obtained from www.flexpart.eu.
- Functionality provided:
- Prepare input 3D-wind fields in hybrid coordinates +
- surface fields for FLEXPART runs
+ ...
"""
# ------------------------------------------------------------------------------
# MODULES
# ------------------------------------------------------------------------------
from gribapi import *
import traceback
-import sys, os
-
+import sys
+import os
# ------------------------------------------------------------------------------
# CLASS
# ------------------------------------------------------------------------------
@@ -54,7 +52,7 @@ class GribTools:
Class for GRIB utilities (new methods) based on GRIB API
'''
# --------------------------------------------------------------------------
- # FUNCTIONS
+ # CLASS FUNCTIONS
# --------------------------------------------------------------------------
def __init__(self, filenames):
'''
diff --git a/python/MARSretrieval.py b/python/MARSretrieval.py
new file mode 100644
index 0000000..e5b6d06
--- /dev/null
+++ b/python/MARSretrieval.py
@@ -0,0 +1,350 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#************************************************************************
+# TODO AP
+#AP
+# -
+# -
+#************************************************************************
+## ------------------------------------------------------------------------------
+# MODULES
+# ------------------------------------------------------------------------------
+import subprocess
+import os
+
+ecapi = True
+try:
+ import ecmwfapi
+except ImportError:
+ ecapi = False
+#from gribapi import *
+# ------------------------------------------------------------------------------
+# CLASS
+# ------------------------------------------------------------------------------
+class MARSretrieval:
+ '''
+ Class for submitting MARS retrievals.
+
+ A description of MARS keywords/arguments and examples of their
+ values can be found here:
+ https://software.ecmwf.int/wiki/display/UDOC/\
+ Identification+keywords#Identificationkeywords-class
+
+ '''
+
+ def __init__(self, server, marsclass = "ei", type = "", levtype = "",
+ levelist = "", repres = "", date = "", resol = "", stream = "",
+ area = "", time = "", step = "", expver = "1", number = "",
+ accuracy = "", grid = "", gaussian = "", target = "",
+ param = ""):
+ '''
+ @Description:
+ Initialises the instance of the MARSretrieval class and
+ defines and assigns a set of the necessary retrieval parameters
+ for the FLEXPART input data.
+ A description of MARS keywords/arguments, their dependencies
+ on each other and examples of their values can be found here:
+
+ https://software.ecmwf.int/wiki/display/UDOC/MARS+keywords
+
+ @Input:
+ self: instance of MARSretrieval
+ For description see class documentation.
+
+ server: instance of ECMWFService (from ECMWF Web-API)
+ This is the connection to the ECMWF data servers.
+ It is needed for the pythonic access of ECMWF data.
+
+ marsclass: string, optional
+ Characterisation of dataset. E.g. EI (ERA-Interim),
+ E4 (ERA40), OD (Operational archive), ea (ERA5).
+ Default is the ERA-Interim dataset "ei".
+
+ type: string, optional
+ Determines the type of fields to be retrieved.
+ Selects between observations, images or fields.
+ Examples for fields: Analysis (an), Forecast (fc),
+ Perturbed Forecast (pf), Control Forecast (cf) and so on.
+ Default is an empty string.
+
+ levtype: string, optional
+ Denotes type of level. Has a direct implication on valid
+ levelist values!
+ E.g. model level (ml), pressure level (pl), surface (sfc),
+ potential vorticity (pv), potential temperature (pt)
+ and depth (dp).
+ Default is an empty string.
+
+ levelist: string, optional
+ Specifies the required levels. It has to have a valid
+ correspondence to the selected levtype.
+ Examples: model level: 1/to/137, pressure levels: 500/to/1000
+ Default is an empty string.
+
+ repres: string, optional
+ Selects the representation of the archived data.
+ E.g. sh - spherical harmonics, gg - Gaussian grid,
+ ll - latitude/longitude, ...
+ Default is an empty string.
+
+ date: string, optional
+ Specifies the Analysis date, the Forecast base date or
+ Observations date. Valid formats are:
+ Absolute as YYYY-MM-DD or YYYYMMDD.
+ Default is an empty string.
+
+ resol: string, optional
+ Specifies the desired triangular truncation of retrieved data,
+ before carrying out any other selected post-processing.
+ The default is automatic truncation (auto), by which the lowest
+ resolution compatible with the value specified in grid is
+ automatically selected for the retrieval.
+ Users wanting to perform post-processing from full spectral
+ resolution should specify Archived Value (av).
+ The following are examples of existing resolutions found in
+ the archive: 63, 106, 159, 213, 255, 319, 399, 511, 799 or 1279.
+ This keyword has no meaning/effect if the archived data is
+ not in spherical harmonics representation.
+ The best selection can be found here:
+ https://software.ecmwf.int/wiki/display/UDOC/\
+ Retrieve#Retrieve-Truncationbeforeinterpolation
+ Default is an empty string.
+
+ stream: string, optional
+ Identifies the forecasting system used to generate the data.
+ E.g. oper (Atmospheric model), enfo (Ensemble forecats), ...
+ Default is an empty string.
+
+ area: string, optional
+ Specifies the desired sub-area of data to be extracted.
+ Areas can be defined to wrap around the globe.
+
+ Latitude values must be given as signed numbers, with:
+ north latitudes (i.e. north of the equator)
+ being positive (e.g: 40.5)
+ south latitutes (i.e. south of the equator)
+ being negative (e.g: -50.5)
+ Longtitude values must be given as signed numbers, with:
+ east longitudes (i.e. east of the 0 degree meridian)
+ being positive (e.g: 35.0)
+ west longitudes (i.e. west of the 0 degree meridian)
+ being negative (e.g: -20.5)
+
+ E.g.: North/West/South/East
+ Default is an empty string.
+
+ time: string, optional
+ Specifies the time of the data in hours and minutes.
+ Valid values depend on the type of data: Analysis time,
+ Forecast base time or First guess verification time
+ (all usually at synoptic hours: 00, 06, 12 and 18 ).
+ Observation time (any combination in hours and minutes is valid,
+ subject to data availability in the archive).
+ The syntax is HHMM or HH:MM. If MM is omitted it defaults to 00.
+ Default is an empty string.
+
+ step: string, optional
+ Specifies the forecast time step from forecast base time.
+ Valid values are hours (HH) from forecast base time. It also
+ specifies the length of the forecast which verifies at
+ First Guess time.
+ E.g. 1/3/6-hourly
+ Default is an empty string.
+
+ expver: string, optional
+ The version of the dataset. Each experiment is assigned a
+ unique code (version). Production data is assigned 1 or 2,
+ and experimental data in Operations 11, 12 ,...
+ Research or Member State's experiments have a four letter
+ experiment identifier.
+ Default is "1".
+
+ number: string, optional
+ Selects the member in ensemble forecast run. (Only then it
+ is necessary.) It has a different meaning depending on
+ the type of data.
+ E.g. Perturbed Forecasts: specifies the Ensemble forecast member
+ Default is an empty string.
+
+ accuracy: string, optional
+ Specifies the number of bits per value to be used in the
+ generated GRIB coded fields.
+ A positive integer may be given to specify the preferred number
+ of bits per packed value. This must not be greater than the
+ number of bits normally used for a Fortran integer on the
+ processor handling the request (typically 32 or 64 bit).
+ Within a compute request the accuracy of the original fields
+ can be passed to the result field by specifying accuracy=av.
+ Default is an empty string.
+
+ grid: string, optional
+ Specifies the output grid which can be either a Gaussian grid
+ or a Latitude/Longitude grid. MARS requests specifying
+ grid=av will return the archived model grid.
+
+ Lat/Lon grid: The grid spacing needs to be an integer
+ fraction of 90 degrees e.g. grid = 0.5/0.5
+
+ Gaussian grid: specified by a letter denoting the type of
+ Gaussian grid followed by an integer (the grid number)
+ representing the number of lines between the Pole and Equator,
+ e.g.
+ grid = F160 - full (or regular) Gaussian grid with
+ 160 latitude lines between the pole and equator
+ grid = N320 - ECMWF original reduced Gaussian grid with
+ 320 latitude lines between the pole and equator,
+ see Reduced Gaussian Grids for grid numbers used at ECMWF
+ grid = O640 - ECMWF octahedral (reduced) Gaussian grid with
+ 640 latitude lines between the pole and equator
+ Default is an empty string.
+
+ gaussian: string, optional
+ This parameter is deprecated and should no longer be used.
+ Specifies the desired type of Gaussian grid for the output.
+ Valid Gaussian grids are quasi-regular (reduced) or regular.
+ Keyword gaussian can only be specified together with
+ keyword grid. Gaussian without grid has no effect.
+ Default is an empty string.
+
+ target: string, optional
+ Specifies a file into which data is to be written after
+ retrieval or manipulation. Path names should always be
+ enclosed in double quotes. The MARS client supports automatic
+ generation of multiple target files using MARS keywords
+ enclosed in square brackets [ ]. If the environment variable
+ MARS_MULTITARGET_STRICT_FORMAT is set to 1 before calling mars,
+ the keyword values will be used in the filename as shown by
+ the ecCodes GRIB tool grib_ls -m, e.g. with
+ MARS_MULTITARGET_STRICT_FORMAT set to 1 the keywords time,
+ expver and param will be formatted as 0600, 0001 and 129.128
+ rather than 600, 1 and 129.
+ Default is an empty string.
+
+ param: string, optional
+ Specifies the meteorological parameter.
+ The list of meteorological parameters in MARS is extensive.
+ Their availability is directly related to their meteorological
+ meaning and, therefore, the rest of directives specified
+ in the MARS request.
+ Meteorological parameters can be specified by their
+ GRIB code (param=130), their mnemonic (param=t) or
+ full name (param=temperature).
+ The list of parameter should be seperated by a "/"-sign.
+ E.g. 130/131/133
+ Default is an empty string.
+
+ @Return:
+ <nothing>
+ '''
+
+ self.server = server
+ self.marsclass = marsclass
+ self.type = type
+ self.levtype = levtype
+ self.levelist = levelist
+ self.repres = repres
+ self.date = date
+ self.resol = resol
+ self.stream = stream
+ self.area = area
+ self.time = time
+ self.step = step
+ self.expver = expver
+ self.number = number
+ self.accuracy = accuracy
+ self.grid = grid
+ self.gaussian = gaussian
+ self.target = target
+ self.param = param
+
+ return
+
+
+ def displayInfo(self):
+ '''
+ @Description:
+ Prints all class attributes and their values.
+
+ @Input:
+ self: instance of MARSretrieval
+ For description see class documentation.
+
+ @Return:
+ <nothing>
+ '''
+ # Get all class attributes and their values as a dictionary
+ attrs = vars(self)
+
+ # iterate through all attributes and print them
+ # with their corresponding values
+ for item in attrs.items():
+ if item[0] in ('server'):
+ pass
+ else:
+ print(item[0] + ': ' + str(item[1]))
+
+ return
+
+ def dataRetrieve(self):
+ '''
+ @Description:
+ Submits a MARS retrieval. Depending on the existence of
+ ECMWF Web-API it is submitted via Python or a
+ subprocess in the Shell. The parameter for the mars retrieval
+ are taken from the defined class attributes.
+
+ @Input:
+ self: instance of MARSretrieval
+ For description see class documentation.
+
+ @Return:
+ <nothing>
+ '''
+ # Get all class attributes and their values as a dictionary
+ attrs = vars(self)
+
+ # convert the dictionary of attributes into a comma
+ # seperated list of attributes with their values
+ # needed for the retrieval call
+ s = 'ret'
+ for k, v in attrs.iteritems():
+ if k in ('server'):
+ continue
+ if k == 'marsclass':
+ k = 'class'
+ if v == '':
+ continue
+ if k.lower() == 'target':
+ target = v
+ else:
+ s = s + ',' + k + '=' + str(v)
+
+ # MARS request via Python script
+ if self.server is not False:
+ try:
+ self.server.execute(s, target)
+ except:
+ print('MARS Request failed, \
+ have you already registered at apps.ecmwf.int?')
+ raise IOError
+ if os.stat(target).st_size == 0:
+ print('MARS Request returned no data - please check request')
+ raise IOError
+ # MARS request via extra process in shell
+ else:
+ s += ',target = "' + target + '"'
+ p = subprocess.Popen(['mars'], stdin=subprocess.PIPE,
+ stdout=subprocess.PIPE,
+ stderr=subprocess.PIPE, bufsize=1)
+ pout = p.communicate(input=s)[0]
+ print(pout.decode())
+
+ if 'Some errors reported' in pout.decode():
+ print('MARS Request failed - please check request')
+ raise IOError
+
+ if os.stat(target).st_size == 0:
+ print('MARS Request returned no data - please check request')
+ raise IOError
+
+ return
diff --git a/python/Tools.py b/python/Tools.py
new file mode 100644
index 0000000..9d55a6b
--- /dev/null
+++ b/python/Tools.py
@@ -0,0 +1,461 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#************************************************************************
+# TODO AP
+#AP
+# -
+#************************************************************************
+"""
+
+"""
+# ------------------------------------------------------------------------------
+# MODULES
+# ------------------------------------------------------------------------------
+from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter
+import os
+import errno
+import sys
+import glob
+from numpy import *
+from gribapi import *
+import Control
+
+# ------------------------------------------------------------------------------
+# FUNCTIONS
+# ------------------------------------------------------------------------------
+
+def interpret_args_and_control():
+ '''
+ @Description:
+ Assigns the command line arguments and reads control file
+ content. Apply default values for non mentioned arguments.
+
+ @Input:
+ <nothing>
+
+ @Return:
+ args: instance of ArgumentParser
+ Contains the commandline arguments from script/program call.
+
+ c: instance of class Control
+ Contains all necessary information of a control file. The parameters
+ are: DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM,
+ NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL, LEVELIST,
+ RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA,
+ SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS,
+ ECFSDIR, MAILOPS, MAILFAIL, GRIB2FLEXPART, DEBUG, INPUTDIR,
+ OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
+ For more information about format and content of the parameter see
+ documentation.
+
+ '''
+ parser = ArgumentParser(description='Retrieve FLEXPART input from \
+ ECMWF MARS archive',
+ formatter_class=ArgumentDefaultsHelpFormatter)
+
+ # the most important arguments
+ parser.add_argument("--start_date", dest="start_date",
+ help="start date YYYYMMDD")
+ parser.add_argument("--end_date", dest="end_date",
+ help="end_date YYYYMMDD")
+ parser.add_argument("--date_chunk", dest="date_chunk", default=None,
+ help="# of days to be retrieved at once")
+
+ # some arguments that override the default in the control file
+ parser.add_argument("--basetime", dest="basetime",
+ help="base such as 00/12 (for half day retrievals)")
+ parser.add_argument("--step", dest="step",
+ help="steps such as 00/to/48")
+ parser.add_argument("--levelist", dest="levelist",
+ help="Vertical levels to be retrieved, e.g. 30/to/60")
+ parser.add_argument("--area", dest="area",
+ help="area defined as north/west/south/east")
+
+ # set the working directories
+ parser.add_argument("--inputdir", dest="inputdir", default=None,
+ help="root directory for storing intermediate files")
+ parser.add_argument("--outputdir", dest="outputdir", default=None,
+ help="root directory for storing output files")
+ parser.add_argument("--flexpart_root_scripts", dest="flexpart_root_scripts",
+ help="FLEXPART root directory (to find grib2flexpart \
+ and COMMAND file)\n\ Normally ECMWFDATA resides in \
+ the scripts directory of the FLEXPART distribution")
+
+ # this is only used by prepareFLEXPART.py to rerun a postprocessing step
+ parser.add_argument("--ppid", dest="ppid",
+ help="Specify parent process id for \
+ rerun of prepareFLEXPART")
+
+ # arguments for job submission to ECMWF, only needed by submit.py
+ parser.add_argument("--job_template", dest='job_template',
+ default="job.temp",
+ help="job template file for submission to ECMWF")
+ parser.add_argument("--queue", dest="queue",
+ help="queue for submission to ECMWF \
+ (e.g. ecgate or cca )")
+ parser.add_argument("--controlfile", dest="controlfile",
+ default='CONTROL.temp',
+ help="file with control parameters")
+ parser.add_argument("--debug", dest="debug", default=0,
+ help="Debug mode - leave temporary files intact")
+
+ args = parser.parse_args()
+
+ # create instance of Control for specified controlfile
+ # and assign the parameters (and default values if necessary)
+ try:
+ c = Control.Control(args.controlfile)
+ except IOError:
+ try:
+ c = Control.Control(localpythonpath + args.controlfile)
+ except:
+ print('Could not read control file "' + args.controlfile + '"')
+ print('Either it does not exist or its syntax is wrong.')
+ print('Try "' + sys.argv[0].split('/')[-1] +
+ ' -h" to print usage information')
+ exit(1)
+
+ # check for having at least a starting date
+ if args.start_date is None and getattr(c, 'start_date') is None:
+ print('start_date specified neither in command line nor \
+ in control file ' + args.controlfile)
+ print('Try "' + sys.argv[0].split('/')[-1] +
+ ' -h" to print usage information')
+ exit(1)
+
+ # save all existing command line parameter to the Control instance
+ # if parameter is not specified through the command line or CONTROL file
+ # set default values
+ if args.start_date is not None:
+ c.start_date = args.start_date
+ if args.end_date is not None:
+ c.end_date = args.end_date
+ if c.end_date is None:
+ c.end_date = c.start_date
+ if args.date_chunk is not None:
+ c.date_chunk = args.date_chunk
+
+ if not hasattr(c, 'debug'):
+ c.debug = args.debug
+
+ if args.inputdir is None and args.outputdir is None:
+ c.inputdir = '../work'
+ c.outputdir = '../work'
+ else:
+ if args.inputdir is not None:
+ c.inputdir = args.inputdir
+ if args.outputdir is None:
+ c.outputdir = args.inputdir
+ if args.outputdir is not None:
+ c.outputdir = args.outputdir
+ if args.inputdir is None:
+ c.inputdir = args.outputdir
+
+ if hasattr(c, 'outputdir') is False and args.outputdir is None:
+ c.outputdir = c.inputdir
+ else:
+ if args.outputdir is not None:
+ c.outputdir = args.outputdir
+
+ if args.area is not None:
+ afloat = '.' in args.area
+ l = args.area.split('/')
+ if afloat:
+ for i in range(len(l)):
+ l[i] = str(int(float(l[i]) * 1000))
+ c.upper, c.left, c.lower, c.right = l
+
+ # NOTE: basetime activates the ''operational mode''
+ if args.basetime is not None:
+ c.basetime = args.basetime
+
+ if args.step is not None:
+ l = args.step.split('/')
+ if 'to' in args.step.lower():
+ if 'by' in args.step.lower():
+ ilist = arange(int(l[0]), int(l[2]) + 1, int(l[4]))
+ c.step = ['{:0>3}'.format(i) for i in ilist]
+ else:
+ myerror(None, args.step + ':\n' +
+ 'please use "by" as well if "to" is used')
+ else:
+ c.step = l
+
+ if args.levelist is not None:
+ c.levelist = args.levelist
+ if 'to' in c.levelist:
+ c.level = c.levelist.split('/')[2]
+ else:
+ c.level = c.levelist.split('/')[-1]
+
+ if args.flexpart_root_scripts is not None:
+ c.flexpart_root_scripts = args.flexpart_root_scripts
+
+ return args, c
+
+
+def cleanup(c):
+ '''
+ @Description:
+ Remove all files from intermediate directory
+ (inputdir from control file).
+
+ @Input:
+ c: instance of class Control
+ Contains all the parameters of control files, which are e.g.:
+ DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
+ STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
+ LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
+ OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
+ ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
+ MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
+ DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
+
+ For more information about format and content of the parameter
+ see documentation.
+
+ @Return:
+ <nothing>
+ '''
+
+ print("cleanup")
+
+ cleanlist = glob.glob(c.inputdir + "/*")
+ for cl in cleanlist:
+ if c.prefix not in cl:
+ silentremove(cl)
+ if c.ecapi is False and (c.ectrans == '1' or c.ecstorage == '1'):
+ silentremove(cl)
+
+ print("Done")
+
+ return
+
+
+def myerror(c, message='ERROR'):
+ '''
+ @Description:
+ Prints a specified error message which can be passed to the function
+ before exiting the program.
+
+ @Input:
+ c: instance of class Control
+ Contains all the parameters of control files, which are e.g.:
+ DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
+ STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
+ LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
+ OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
+ ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
+ MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
+ DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
+
+ For more information about format and content of the parameter
+ see documentation.
+
+ message: string, optional
+ Error message. Default value is "ERROR".
+
+ @Return:
+ <nothing>
+ '''
+ # uncomment if user wants email notification directly from python
+ #try:
+ #target = c.mailfail
+ #except AttributeError:
+ #target = os.getenv('USER')
+
+ #if(type(target) is not list):
+ #target = [target]
+
+ print(message)
+
+ # uncomment if user wants email notification directly from python
+ #for t in target:
+ #p = subprocess.Popen(['mail','-s ECMWFDATA v7.0 ERROR', os.path.expandvars(t)],
+ # stdin = subprocess.PIPE, stdout = subprocess.PIPE,
+ # stderr = subprocess.PIPE, bufsize = 1)
+ #tr = '\n'.join(traceback.format_stack())
+ #pout = p.communicate(input = message+'\n\n'+tr)[0]
+ #print 'Email sent to '+os.path.expandvars(t) # +' '+pout.decode()
+
+ exit(1)
+
+ return
+
+
+def normalexit(c, message='Done!'):
+ '''
+ @Description:
+ Prints a specific exit message which can be passed to the function.
+
+ @Input:
+ c: instance of class Control
+ Contains all the parameters of control files, which are e.g.:
+ DAY1(start_date), DAY2(end_date), DTIME, MAXSTEP, TYPE, TIME,
+ STEP, CLASS(marsclass), STREAM, NUMBER, EXPVER, GRID, LEFT,
+ LOWER, UPPER, RIGHT, LEVEL, LEVELIST, RESOL, GAUSS, ACCURACY,
+ OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, SMOOTH, FORMAT,
+ ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR,
+ MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR, BASETIME
+ DATE_CHUNK, DEBUG, INPUTDIR, OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
+
+ For more information about format and content of the parameter
+ see documentation.
+
+ message: string, optional
+ Message for exiting program. Default value is "Done!".
+
+ @Return:
+ <nothing>
+
+ '''
+ # Uncomment if user wants notification directly from python
+ #try:
+ #target = c.mailops
+ #if(type(target) is not list):
+ #target = [target]
+ #for t in target:
+ #p = subprocess.Popen(['mail','-s ECMWFDATA v7.0 normal exit',
+ # os.path.expandvars(t)],
+ # stdin = subprocess.PIPE,
+ # stdout = subprocess.PIPE,
+ # stderr = subprocess.PIPE, bufsize = 1)
+ #pout = p.communicate(input = message+'\n\n')[0]
+ #print pout.decode()
+ #except:
+ #pass
+
+ print(message)
+
+ return
+
+
+def product(*args, **kwds):
+ '''
+ @Description:
+ This method is taken from an example at the ECMWF wiki website.
+ https://software.ecmwf.int/wiki/display/GRIB/index.py; 2018-03-16
+
+ This method combines the single characters of the passed arguments
+ with each other. So that each character of each argument value
+ will be combined with each character of the other arguments as a tuple.
+
+ Example:
+ product('ABCD', 'xy') --> Ax Ay Bx By Cx Cy Dx Dy
+ product(range(2), repeat = 3) --> 000 001 010 011 100 101 110 111
+
+ @Input:
+ *args: tuple
+ Positional arguments (arbitrary number).
+
+ **kwds: dictionary
+ Contains all the keyword arguments from *args.
+
+ @Return:
+ prod: tuple
+ Return will be done with "yield". A tuple of combined arguments.
+ See example in description above.
+ '''
+
+ pools = map(tuple, args) * kwds.get('repeat', 1)
+ result = [[]]
+ for pool in pools:
+ result = [x + [y] for x in result for y in pool]
+ for prod in result:
+ yield tuple(prod)
+
+ return
+
+
+def silentremove(filename):
+ '''
+ @Description:
+ Removes the file which name is passed to the function if
+ it exists. The function does not fail if the file does not
+ exist.
+
+ @Input:
+ filename: string
+ The name of the file to be removed without notification.
+
+ @Return:
+ <nothing>
+ '''
+ try:
+ os.remove(filename)
+ except OSError as e:
+ # this would be "except OSError, e:" before Python 2.6
+ if e.errno is not errno.ENOENT:
+ # errno.ENOENT = no such file or directory
+ raise # re-raise exception if a different error occured
+
+ return
+
+
+def init128(fn):
+ '''
+ @Description:
+ Opens and reads the grib file with table 128 information.
+
+ @Input:
+ fn: string
+ Path to file of ECMWF grib table number 128.
+
+ @Return:
+ table128: dictionary
+ Contains the ECMWF grib table 128 information.
+ The key is the parameter number and the value is the
+ short name of the parameter.
+ '''
+ table128 = dict()
+ with open(fn) as f:
+ fdata = f.read().split('\n')
+ for data in fdata:
+ if data[0] != '!':
+ table128[data[0:3]] = data[59:64].strip()
+
+ return table128
+
+
+def toparamId(pars, table):
+ '''
+ @Description:
+ Transform parameter names to parameter ids
+ with ECMWF grib table 128.
+
+ @Input:
+ pars: string
+ Addpar argument from control file in the format of
+ parameter names instead of ids. The parameter short
+ names are sepearted with "/" and they are passed as
+ one single string.
+
+ table: dictionary
+ Contains the ECMWF grib table 128 information.
+ The key is the parameter number and the value is the
+ short name of the parameter.
+
+ @Return:
+ ipar: list of integer
+ List of addpar parameters from control file transformed to
+ parameter ids in the format of integer.
+ '''
+ cpar = pars.upper().split('/')
+ ipar = []
+ for par in cpar:
+ found = False
+ for k, v in table.iteritems():
+ if par == k or par == v:
+ ipar.append(int(k))
+ found = True
+ break
+ if found is False:
+ print('Warning: par ' + par + ' not found in table 128')
+
+ return ipar
+
+def getListAsString(listobj):
+ '''
+ @Description:
+ '''
+ return ", ".join( str(l) for l in listobj)
\ No newline at end of file
diff --git a/python/UIOTools.py b/python/UIOFiles.py
similarity index 59%
rename from python/UIOTools.py
rename to python/UIOFiles.py
index 3ee9046..faeead6 100644
--- a/python/UIOTools.py
+++ b/python/UIOFiles.py
@@ -3,10 +3,9 @@
#************************************************************************
# TODO AP
#AP
-# - File name und Klassenname gleichsetzen?
# - checken welche regelmässigen methoden auf diese Files noch angewendet werden
# und dann hier implementieren
-# - löschen?
+# - add description of file!
#************************************************************************
"""
@Author: Anne Fouilloux (University of Oslo)
@@ -14,39 +13,39 @@
@Date: October 2014
@ChangeHistory:
- February 2018 - Anne Philipp (University of Vienna):
+ November 2015 - Leopold Haimberger (University of Vienna):
+ - modified method listFiles to work with glob instead of listdir
+ - added pattern search in method listFiles
+
+ February 2018 - Anne Philipp (University of Vienna):
- applied PEP8 style guide
- added documentation
+ - optimisation of method listFiles since it didn't work correctly
+ for sub directories
+ - additional speed up of method listFiles
@License:
- (C) Copyright 2014 UIO.
+ (C) Copyright 2014-2018.
This software is licensed under the terms of the Apache Licence Version 2.0
which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
@Requirements:
- - A standard python 2.6 or 2.7 installation
- - dateutils
- - matplotlib (optional, for debugging)
- - ECMWF specific packages, all available from https://software.ecmwf.int/
- ECMWF WebMARS, gribAPI with python enabled, emoslib and
- ecaccess web toolkit
+ A standard python 2.6 or 2.7 installation
@Description:
- Further documentation may be obtained from www.flexpart.eu.
-
- Functionality provided:
- Prepare input 3D-wind fields in hybrid coordinates +
- surface fields for FLEXPART runs
+ ...
"""
# ------------------------------------------------------------------------------
# MODULES
# ------------------------------------------------------------------------------
import os
import glob
-
+import fnmatch
+import time
+import profiling
# ------------------------------------------------------------------------------
-# Class
+# CLASS
# ------------------------------------------------------------------------------
class UIOFiles:
'''
@@ -55,7 +54,7 @@ class UIOFiles:
with the instance of the class.
'''
# --------------------------------------------------------------------------
- # FUNCTIONS
+ # CLASS FUNCTIONS
# --------------------------------------------------------------------------
def __init__(self, suffix):
'''
@@ -79,7 +78,8 @@ class UIOFiles:
return
- def listFiles(self, path, pattern):
+ #@profiling.timefn
+ def listFiles(self, path, pattern, callid=0):
'''
@Description:
Lists all files in the directory with the matching
@@ -97,37 +97,35 @@ class UIOFiles:
Regular expression pattern. For example:
'*OG_acc_SL*.'+c.ppid+'.*'
+ callid: integer
+ Id which tells the function if its the first call
+ or a recursive call. Default and first call is 0.
+ Everything different from 0 is ment to be a recursive case.
+
@Return:
<nothing>
'''
+ # initialize variable in first function call
+ if callid == 0:
+ self.files = []
+
# Get the absolute path
path = os.path.abspath(path)
- # Get a list of files in pathname
- filesInCurDir0 = glob.glob(path + '/' + pattern)
- filesInCurDir = []
- for f in filesInCurDir0:
- filesInCurDir.append(f.split('/')[-1])
-
- self.counter = 0
- self.files = []
- # Traverse through all files
- for file in filesInCurDir:
- curFile = os.path.join(path, file)
-
- # Check if it's a normal file or directory
- if os.path.isfile(curFile):
- # Get the file extension
- fileNoExt, curFileExtension = os.path.splitext(curFile)
- # Check if the file has an extension of typical video files
- if curFileExtension in self.suffix:
- # We have got a file file! Increment the counter
- self.counter += 1
- # add this filename in the list
- self.files.append(curFile)
- else:
- # We got a directory, enter into it for further processing
- self.listFiles(curFile)
-
- return
\ No newline at end of file
+ # get the file list of the path if its not a directory and
+ # if it contains one of the suffixes
+ self.files.extend([os.path.join(path, k) for k in os.listdir(path)
+ if fnmatch.fnmatch(k, pattern) and
+ os.path.splitext(k)[-1] in self.suffix])
+
+ # find possible sub-directories in the path
+ subdirs = [s for s in os.listdir(path)
+ if os.path.isdir(os.path.join(path, s))]
+
+ # do recursive calls for sub-direcorties
+ if subdirs:
+ for subdir in subdirs:
+ self.listFiles(os.path.join(path, subdir), pattern, callid=1)
+
+ return
diff --git a/python/compilejob.ksh b/python/compilejob.ksh
index c943665..d276521 100644
--- a/python/compilejob.ksh
+++ b/python/compilejob.ksh
@@ -4,7 +4,7 @@
# start with ecaccess-job-submit -queueName ecgb NAME_OF_THIS_FILE on gateway server
# start with sbatch NAME_OF_THIS_FILE directly on machine
-#SBATCH --workdir=/scratch/ms/spatlh00/lh0
+#SBATCH --workdir=/scratch/ms/at/km4a
#SBATCH --qos=normal
#SBATCH --job-name=flex_ecmwf
#SBATCH --output=flex_ecmwf.%j.out
@@ -24,7 +24,7 @@
##PBS -l EC_memory_per_task=3200MB
set -x
-export VERSION=7.0
+export VERSION=7.1
case $HOST in
*ecg*)
module load python
@@ -32,8 +32,8 @@ case $HOST in
module unload emos
module load grib_api/1.14.5
module load emos/437-r64
-export FLEXPART_ROOT_SCRIPTS=$HOME
-# export ECMWFDATA=$FLEXPART_ROOT/ECMWFDATA7.0
+export FLEXPART_ROOT_SCRIPTS=${HOME}
+# export ECMWFDATA=$FLEXPART_ROOT/ECMWFDATA$VERSION
# export PYTHONPATH=$ECMWFDATA/python
# export PATH=${PATH}:$ECMWFDATA/python
export MAKEFILE=Makefile.gfortran
@@ -48,8 +48,8 @@ export FLEXPART_ROOT_SCRIPTS=$HOME
echo $HOME | awk -F / '{print $1, $2, $3, $4}'
export GROUP=`echo $HOME | awk -F / '{print $4}'`
export SCRATCH=/scratch/ms/${GROUP}/${USER}
-export FLEXPART_ROOT_SCRIPTS=$HOME
-# export ECMWFDATA=$FLEXPART_ROOT/ECMWFDATA7.0
+export FLEXPART_ROOT_SCRIPTS=${HOME}
+# export ECMWFDATA=$FLEXPART_ROOT/ECMWFDATA$VERSION
# export PYTHONPATH=$ECMWFDATA/python
# export PATH=${PATH}:$ECMWFDATA/python
export MAKEFILE=Makefile.CRAY
@@ -58,7 +58,7 @@ esac
mkdir -p $FLEXPART_ROOT_SCRIPTS/ECMWFDATA$VERSION
cd $FLEXPART_ROOT_SCRIPTS/ECMWFDATA$VERSION # if FLEXPART_ROOT is not set this means cd to the home directory
-tar -xvf $SCRATCH/ECMWFDATA$VERSION.tar
+tar -xvf $HOME/ECMWFDATA$VERSION.tar
cd src
\rm *.o *.mod CONVERT2
make -f $MAKEFILE >flexcompile 2>flexcompile
diff --git a/python/compilejob.temp b/python/compilejob.temp
index d2da195..e699a01 100644
--- a/python/compilejob.temp
+++ b/python/compilejob.temp
@@ -24,7 +24,7 @@
##PBS -l EC_memory_per_task=3200MB
set -x
-export VERSION=7.0
+export VERSION=7.1
case $HOST in
*ecg*)
module load python
@@ -33,7 +33,7 @@ case $HOST in
module load grib_api/1.14.5
module load emos/437-r64
export FLEXPART_ROOT_SCRIPTS=
-# export ECMWFDATA=$FLEXPART_ROOT/ECMWFDATA7.0
+# export ECMWFDATA=$FLEXPART_ROOT/ECMWFDATA$VERSION
# export PYTHONPATH=$ECMWFDATA/python
# export PATH=${PATH}:$ECMWFDATA/python
export MAKEFILE=Makefile.gfortran
@@ -49,7 +49,7 @@ case $HOST in
export GROUP=`echo $HOME | awk -F / '{print $4}'`
export SCRATCH=/scratch/ms/${GROUP}/${USER}
export FLEXPART_ROOT_SCRIPTS=
-# export ECMWFDATA=$FLEXPART_ROOT/ECMWFDATA7.0
+# export ECMWFDATA=$FLEXPART_ROOT/ECMWFDATA$VERSION
# export PYTHONPATH=$ECMWFDATA/python
# export PATH=${PATH}:$ECMWFDATA/python
export MAKEFILE=Makefile.CRAY
@@ -58,7 +58,7 @@ esac
mkdir -p $FLEXPART_ROOT_SCRIPTS/ECMWFDATA$VERSION
cd $FLEXPART_ROOT_SCRIPTS/ECMWFDATA$VERSION # if FLEXPART_ROOT is not set this means cd to the home directory
-tar -xvf $SCRATCH/ECMWFDATA$VERSION.tar
+tar -xvf $HOME/ECMWFDATA$VERSION.tar
cd src
\rm *.o *.mod CONVERT2
make -f $MAKEFILE >flexcompile 2>flexcompile
diff --git a/python/getMARSdata.py b/python/getMARSdata.py
old mode 100644
new mode 100755
index e319a8f..278176d
--- a/python/getMARSdata.py
+++ b/python/getMARSdata.py
@@ -1,36 +1,51 @@
#!/usr/bin/env python
-#
-# This software is licensed under the terms of the Apache Licence Version 2.0
-# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
-#
-# Functionality provided: Prepare input 3D-wind fields in hybrid coordinates + surface fields for FLEXPART runs
-#
-# Creation: October 2014 - Anne Fouilloux - University of Oslo
-# Extension November 2015 - Leopold Haimberger - University of Vienna for:
-# - using the WebAPI also for general MARS retrievals
-# - job submission on ecgate and cca
-# - job templates suitable for twice daily operational dissemination
-# - dividing retrievals of longer periods into digestable chunks
-# - retrieve also longer term forecasts, not only analyses and short term forecast data
-# - conversion into GRIB2
-# - conversion into .fp format for faster execution of FLEXPART
-#
-#
-# Further documentation may be obtained from www.flexpart.eu
-#
-# Requirements:
-# in addition to a standard python 2.6 or 2.7 installation the following packages need to be installed
-# ECMWF WebMARS, gribAPI with python enabled, emoslib, ecaccess web toolkit, all available from https://software.ecmwf.int/
-# dateutils
-# matplotlib (optional, for debugging)
-#
-# Get MARS GRIB fields from ECMWF for FLEXPART
-#
-
-#import socket
-
-#hostname=socket.gethostname()
-#ecapi= 'ecmwf' not in hostname
+# -*- coding: utf-8 -*-
+#************************************************************************
+# TODO AP
+# - Change History ist nicht angepasst ans File!
+# - add file description
+#************************************************************************
+"""
+@Author: Anne Fouilloux (University of Oslo)
+
+@Date: October 2014
+
+@ChangeHistory:
+ November 2015 - Leopold Haimberger (University of Vienna):
+ - using the WebAPI also for general MARS retrievals
+ - job submission on ecgate and cca
+ - job templates suitable for twice daily operational dissemination
+ - dividing retrievals of longer periods into digestable chunks
+ - retrieve also longer term forecasts, not only analyses and
+ short term forecast data
+ - conversion into GRIB2
+ - conversion into .fp format for faster execution of FLEXPART
+
+ February 2018 - Anne Philipp (University of Vienna):
+ - applied PEP8 style guide
+ - added documentation
+ - minor changes in programming style for consistence
+
+@License:
+ (C) Copyright 2014-2018.
+
+ This software is licensed under the terms of the Apache Licence Version 2.0
+ which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+
+@Requirements:
+ - A standard python 2.6 or 2.7 installation
+ - dateutils
+ - ECMWF specific packages, all available from https://software.ecmwf.int/
+ ECMWF WebMARS, gribAPI with python enabled, emoslib and
+ ecaccess web toolkit
+
+@Description:
+ Further documentation may be obtained from www.flexpart.eu.
+
+"""
+# ------------------------------------------------------------------------------
+# MODULES
+# ------------------------------------------------------------------------------
try:
ecapi=True
import ecmwfapi
@@ -41,19 +56,23 @@ import calendar
import shutil
import datetime
import time
-import os,glob,sys,inspect
-#from string import strip
-from argparse import ArgumentParser,ArgumentDefaultsHelpFormatter
+import os
+import glob
+import sys
+import inspect
# add path to submit.py to pythonpath so that python finds its buddies
localpythonpath=os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
if localpythonpath not in sys.path:
sys.path.append(localpythonpath)
-from FlexpartTools import ECFlexpart, \
- Control, myerror, normalexit, \
+from Control import Control
+from Tools import myerror, normalexit, \
interpret_args_and_control
+from ECFlexpart import ECFlexpart
-
+# ------------------------------------------------------------------------------
+# FUNCTION
+# ------------------------------------------------------------------------------
def getMARSdata(args, c):
@@ -70,7 +89,7 @@ def getMARSdata(args, c):
c.ecapi = ecapi
print 'ecapi:', c.ecapi
-# Retrieve ERA interim data for running flexpart
+# Retrieve EC data for running flexpart
#AP change this variant to correct format conversion with datetime
#AP import datetime and timedelta explicitly
syear = int(c.start_date[:4])
@@ -80,6 +99,7 @@ def getMARSdata(args, c):
startm1 = start - datetime.timedelta(days=1)
if c.basetime == '00':
start = startm1
+
eyear = int(c.end_date[:4])
emonth = int(c.end_date[4:6])
eday = int(c.end_date[6:])
@@ -90,66 +110,103 @@ def getMARSdata(args, c):
endp1 = end + datetime.timedelta(days=2)
datechunk = datetime.timedelta(days=int(c.date_chunk))
+
+ # retrieving of accumulated data fields (flux data), (maximum one month)
+
+ # remove old files
print 'removing content of ' + c.inputdir
- tobecleaned = glob.glob(c.inputdir + '/*_acc_*.' + str(os.getppid()) + '.*.grb')
+ tobecleaned = glob.glob(c.inputdir + '/*_acc_*.' + \
+ str(os.getppid()) + '.*.grb')
for f in tobecleaned:
os.remove(f)
- times=None
- if c.maxstep<24:
- day=startm1
- while day<endp1:
- # we need to retrieve MARS data for this period (maximum one month)
- flexpart = ECFlexpart(c,fluxes=True)
- if day+datechunk-datetime.timedelta(days=1)<endp1:
- dates= day.strftime("%Y%m%d") + "/to/" + (day+datechunk-datetime.timedelta(days=1)).strftime("%Y%m%d")
+ times = None
+ # if forecast for maximum one day (upto 23h) are to be retrieved,
+ # collect accumulation data (flux data)
+ # with additional days in the beginning and at the end
+ # (used for complete disaggregation of original period)
+ if c.maxstep < 24:
+ day = startm1
+ while day < endp1:
+ # retrieve MARS data for the whole period
+ flexpart = ECFlexpart(c, fluxes=True)
+ tmpday = day + datechunk - datetime.timedelta(days=1)
+ if tmpday < endp1:
+ dates = day.strftime("%Y%m%d") + "/to/" + \
+ tmpday.strftime("%Y%m%d")
else:
- dates= day.strftime("%Y%m%d") + "/to/" + end.strftime("%Y%m%d")
+ dates = day.strftime("%Y%m%d") + "/to/" + \
+ end.strftime("%Y%m%d")
print "retrieve " + dates + " in dir " + c.inputdir
+
try:
flexpart.retrieve(server, dates, times, c.inputdir)
except IOError:
myerror(c,'MARS request failed')
- day+=datechunk
+ day += datechunk
+
+ # if forecast data longer than 24h are to be retrieved,
+ # collect accumulation data (flux data)
+ # with the exact start and end date
+ # (disaggregation will be done for the
+ # exact time period with boundary conditions)
else:
- day=start
- while day<=end:
- # we need to retrieve MARS data for this period (maximum one month)
- flexpart = ECFlexpart(c,fluxes=True)
- if day+datechunk-datetime.timedelta(days=1)<end:
- dates= day.strftime("%Y%m%d") + "/to/" + (day+datechunk-datetime.timedelta(days=1)).strftime("%Y%m%d")
+ day = start
+ while day <= end:
+ # retrieve MARS data for the whole period
+ flexpart = ECFlexpart(c, fluxes=True)
+ tmpday = day + datechunk - datetime.timedelta(days=1)
+ if tmpday < end:
+ dates = day.strftime("%Y%m%d") + "/to/" + \
+ tmpday.trftime("%Y%m%d")
else:
- dates= day.strftime("%Y%m%d") + "/to/" + end.strftime("%Y%m%d")
+ dates = day.strftime("%Y%m%d") + "/to/" + \
+ end.strftime("%Y%m%d")
print "retrieve " + dates + " in dir " + c.inputdir
- flexpart.retrieve(server, dates, times, c.inputdir)
- day+=datechunk
+ try:
+ flexpart.retrieve(server, dates, times, c.inputdir)
+ except IOError:
+ myerror(c, 'MARS request failed')
+
+ day += datechunk
+ # retrieving of normal data fields (non flux data), (maximum one month)
- tobecleaned=glob.glob(c.inputdir+'/*__*.'+str(os.getppid())+'.*.grb')
+ # remove old *__* files
+ tobecleaned = glob.glob(c.inputdir + '/*__*.' +
+ str(os.getppid()) + '.*.grb')
for f in tobecleaned:
os.remove(f)
- day=start
- times=None
- while day<=end:
-
- # we need to retrieve MARS data for this period (maximum one month)
- flexpart = ECFlexpart(c)
- if day+datechunk-datetime.timedelta(days=1)<end:
- dates= day.strftime("%Y%m%d") + "/to/" + (day+datechunk-datetime.timedelta(days=1)).strftime("%Y%m%d")
+ day = start
+ times = None
+ while day <= end:
+ # retrieve MARS data for the whole period
+ flexpart = ECFlexpart(c, fluxes=False)
+ tmpday = day+datechunk-datetime.timedelta(days=1)
+ if tmpday < end:
+ dates = day.strftime("%Y%m%d") + "/to/" + \
+ tmpday.strftime("%Y%m%d")
else:
- dates= day.strftime("%Y%m%d") + "/to/" + end.strftime("%Y%m%d")
+ dates = day.strftime("%Y%m%d") + "/to/" + \
+ end.strftime("%Y%m%d")
+
print "retrieve " + dates + " in dir " + c.inputdir
- flexpart.retrieve(server, dates, times, c.inputdir)
- day+=datechunk
+ try:
+ flexpart.retrieve(server, dates, times, c.inputdir)
+ except IOError:
+ myerror(c, 'MARS request failed')
+
+ day += datechunk
+ return
if __name__ == "__main__":
- args,c=interpret_args_and_control()
- getMARSdata(args,c)
+ args, c = interpret_args_and_control()
+ getMARSdata(args, c)
normalexit(c)
diff --git a/python/install.py b/python/install.py
old mode 100644
new mode 100755
index 7385f9d..d73b08c
--- a/python/install.py
+++ b/python/install.py
@@ -63,17 +63,124 @@ import inspect
# add path to submit.py to pythonpath so that python finds its buddies
localpythonpath=os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
sys.path.append(localpythonpath)
-from UIOTools import UIOFiles
+from UIOFiles import UIOFiles
from string import strip
from argparse import ArgumentParser,ArgumentDefaultsHelpFormatter
from GribTools import GribTools
-from FlexpartTools import ECFlexpart, Control, install_args_and_control
+from Control import Control
from getMARSdata import getMARSdata
from prepareFLEXPART import prepareFLEXPART
+from ECFlexpart import ECFlexpart
# ------------------------------------------------------------------------------
# FUNCTIONS
# ------------------------------------------------------------------------------
+def install_args_and_control():
+ '''
+ @Description:
+ Assigns the command line arguments for installation and reads
+ control file content. Apply default values for non mentioned arguments.
+
+ @Input:
+ <nothing>
+
+ @Return:
+ args: instance of ArgumentParser
+ Contains the commandline arguments from script/program call.
+
+ c: instance of class Control
+ Contains all necessary information of a control file. The parameters
+ are: DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM,
+ NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL, LEVELIST,
+ RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA,
+ SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS,
+ ECFSDIR, MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR
+ For more information about format and content of the parameter see
+ documentation.
+ '''
+ parser = ArgumentParser(description='Install ECMWFDATA software locally or \
+ on ECMWF machines',
+ formatter_class=ArgumentDefaultsHelpFormatter)
+
+ parser.add_argument('--target', dest='install_target',
+ help="Valid targets: local | ecgate | cca , \
+ the latter two are at ECMWF")
+ parser.add_argument("--makefile", dest="makefile",
+ help='Name of Makefile to use for compiling CONVERT2')
+ parser.add_argument("--ecuid", dest="ecuid",
+ help='user id at ECMWF')
+ parser.add_argument("--ecgid", dest="ecgid",
+ help='group id at ECMWF')
+ parser.add_argument("--gateway", dest="gateway",
+ help='name of local gateway server')
+ parser.add_argument("--destination", dest="destination",
+ help='ecaccess destination, e.g. leo@genericSftp')
+
+ parser.add_argument("--flexpart_root_scripts", dest="flexpart_root_scripts",
+ help="FLEXPART root directory on ECMWF servers \
+ (to find grib2flexpart and COMMAND file)\n\
+ Normally ECMWFDATA resides in the scripts directory \
+ of the FLEXPART distribution, thus the:")
+
+# arguments for job submission to ECMWF, only needed by submit.py
+ parser.add_argument("--job_template", dest='job_template',
+ default="job.temp.o",
+ help="job template file for submission to ECMWF")
+
+ parser.add_argument("--controlfile", dest="controlfile",
+ default='CONTROL.temp',
+ help="file with control parameters")
+
+ args = parser.parse_args()
+
+ try:
+ c = Control(args.controlfile)
+ except:
+ print('Could not read control file "' + args.controlfile + '"')
+ print('Either it does not exist or its syntax is wrong.')
+ print('Try "' + sys.argv[0].split('/')[-1] +
+ ' -h" to print usage information')
+ exit(1)
+
+ if args.install_target != 'local':
+ if (args.ecgid is None or args.ecuid is None or args.gateway is None
+ or args.destination is None):
+ print('Please enter your ECMWF user id and group id as well as \
+ the \nname of the local gateway and the ectrans \
+ destination ')
+ print('with command line options --ecuid --ecgid \
+ --gateway --destination')
+ print('Try "' + sys.argv[0].split('/')[-1] +
+ ' -h" to print usage information')
+ print('Please consult ecaccess documentation or ECMWF user support \
+ for further details')
+ sys.exit(1)
+ else:
+ c.ecuid = args.ecuid
+ c.ecgid = args.ecgid
+ c.gateway = args.gateway
+ c.destination = args.destination
+
+ try:
+ c.makefile = args.makefile
+ except:
+ pass
+
+ if args.install_target == 'local':
+ if args.flexpart_root_scripts is None:
+ c.flexpart_root_scripts = '../'
+ else:
+ c.flexpart_root_scripts = args.flexpart_root_scripts
+
+ if args.install_target != 'local':
+ if args.flexpart_root_scripts is None:
+ c.ec_flexpart_root_scripts = '${HOME}'
+ else:
+ c.ec_flexpart_root_scripts = args.flexpart_root_scripts
+
+ return args, c
+
+
def main():
'''
'''
@@ -138,7 +245,7 @@ def install_via_gateway(c, target):
data = '##PBS -o /scratch/ms/' + c.ecgid + '/' + \
c.ecuid + 'flex_ecmwf.$Jobname.$Job_ID.out'
if 'export PATH=${PATH}:' in data:
- data += c.ec_flexpart_root_scripts + '/ECMWFDATA7.0/python'
+ data += c.ec_flexpart_root_scripts + '/ECMWFDATA7.1/python'
if 'cat>>' in data or 'cat >>' in data:
i = data.index('>')
fo.write(data[:i] + data[i+1:] + '\n')
@@ -170,21 +277,21 @@ def install_via_gateway(c, target):
if os.path.abspath(ecd) != os.path.abspath(c.flexpart_root_scripts):
os.chdir('/')
p = subprocess.check_call(['tar', '-cvf',
- ecd + '../ECMWFDATA7.0.tar',
+ ecd + '../ECMWFDATA7.1.tar',
ecd + 'python',
ecd + 'grib_templates',
ecd + 'src'])
try:
- os.makedirs(c.flexpart_root_scripts + '/ECMWFDATA7.0')
+ os.makedirs(c.flexpart_root_scripts + '/ECMWFDATA7.1')
except:
pass
- os.chdir(c.flexpart_root_scripts + '/ECMWFDATA7.0')
+ os.chdir(c.flexpart_root_scripts + '/ECMWFDATA7.1')
p = subprocess.check_call(['tar', '-xvf',
- ecd + '../ECMWFDATA7.0.tar'])
- os.chdir(c.flexpart_root_scripts + '/ECMWFDATA7.0/src')
+ ecd + '../ECMWFDATA7.1.tar'])
+ os.chdir(c.flexpart_root_scripts + '/ECMWFDATA7.1/src')
os.chdir('../src')
- print(('install ECMWFDATA7.0 software on ' + target + ' in directory '
+ print(('install ECMWFDATA7.1 software on ' + target + ' in directory '
+ os.getcwd()))
if c.makefile is None:
makefile = 'Makefile.local.ifort'
@@ -196,27 +303,27 @@ def install_via_gateway(c, target):
try:
print(('Using makefile: ' + makefile))
p = subprocess.check_call(['make', '-f', makefile])
- p = subprocess.check_call(['ls', '-l',' CONVERT2'])
+ p = subprocess.check_call(['ls', '-l','CONVERT2'])
except:
print('compile failed - please edit ' + makefile +
' or try another Makefile in the src directory.')
- print('most likely GRIB_API_INCLUDE_DIR, GRIB_API_LIB \
- and EMOSLIB must be adapted.')
+ print('most likely GRIB_API_INCLUDE_DIR, GRIB_API_LIB '
+ 'and EMOSLIB must be adapted.')
print('Available Makefiles:')
print(glob.glob('Makefile*'))
elif target.lower() == 'ecgate':
os.chdir('/')
p = subprocess.check_call(['tar', '-cvf',
- ecd + '../ECMWFDATA7.0.tar',
+ ecd + '../ECMWFDATA7.1.tar',
ecd + 'python',
ecd + 'grib_templates',
ecd + 'src'])
try:
p = subprocess.check_call(['ecaccess-file-put',
- ecd + '../ECMWFDATA7.0.tar',
- 'ecgate:/scratch/ms/' + c.ecgid + '/' +
- c.ecuid + '/ECMWFDATA7.0.tar'])
+ ecd + '../ECMWFDATA7.1.tar',
+ 'ecgate:/home/ms/' + c.ecgid + '/' +
+ c.ecuid + '/ECMWFDATA7.1.tar'])
except:
print('ecaccess-file-put failed! Probably the eccert key has expired.')
exit(1)
@@ -224,37 +331,37 @@ def install_via_gateway(c, target):
'-queueName',
target,
ecd + 'python/compilejob.ksh'])
- print('compilejob.ksh has been submitted to ecgate for \
- installation in ' + c.ec_flexpart_root_scripts +
- '/ECMWFDATA7.0')
- print('You should get an email with subject flexcompile within \
- the next few minutes')
+ print('compilejob.ksh has been submitted to ecgate for '
+ 'installation in ' + c.ec_flexpart_root_scripts +
+ '/ECMWFDATA7.1')
+ print('You should get an email with subject flexcompile within '
+ 'the next few minutes')
elif target.lower() == 'cca':
os.chdir('/')
p = subprocess.check_call(['tar', '-cvf',
- ecd + '../ECMWFDATA7.0.tar',
+ ecd + '../ECMWFDATA7.1.tar',
ecd + 'python',
ecd + 'grib_templates',
ecd + 'src'])
try:
p = subprocess.check_call(['ecaccess-file-put',
- ecd + '../ECMWFDATA7.0.tar',
- 'cca:/scratch/ms/' + c.ecgid + '/' +
- c.ecuid + '/ECMWFDATA7.0.tar'])
+ ecd + '../ECMWFDATA7.1.tar',
+ 'cca:/home/ms/' + c.ecgid + '/' +
+ c.ecuid + '/ECMWFDATA7.1.tar'])
except:
- print('ecaccess-file-put failed! \
- Probably the eccert key has expired.')
+ print('ecaccess-file-put failed! '
+ 'Probably the eccert key has expired.')
exit(1)
p=subprocess.check_call(['ecaccess-job-submit',
'-queueName',
target,
- ecd + 'python/compilejob.ksh']))
+ ecd + 'python/compilejob.ksh'])
print('compilejob.ksh has been submitted to cca for installation in ' +
- c.ec_flexpart_root_scripts + '/ECMWFDATA7.0')
- print('You should get an email with subject flexcompile \
- within the next few minutes')
+ c.ec_flexpart_root_scripts + '/ECMWFDATA7.1')
+ print('You should get an email with subject flexcompile '
+ 'within the next few minutes')
else:
print('ERROR: unknown installation target ', target)
diff --git a/python/job.ksh b/python/job.ksh
index acc7256..70e561c 100644
--- a/python/job.ksh
+++ b/python/job.ksh
@@ -4,7 +4,7 @@
# start with ecaccess-job-submit -queueName ecgb NAME_OF_THIS_FILE on gateway server
# start with sbatch NAME_OF_THIS_FILE directly on machine
-#SBATCH --workdir=/scratch/ms/spatlh00/lh0
+#SBATCH --workdir=/scratch/ms/at/km4a
#SBATCH --qos=normal
#SBATCH --job-name=flex_ecmwf
#SBATCH --output=flex_ecmwf.%j.out
@@ -24,7 +24,7 @@
##PBS -l EC_memory_per_task=32000MB
set -x
-
+export VERSION=7.1
case $HOST in
*ecg*)
module load python
@@ -32,9 +32,9 @@ case $HOST in
module unload emos
module load grib_api/1.14.5
module load emos/437-r64
-# export ECMWFDATA=$HOME/ECMWFDATA7.0
+# export ECMWFDATA=$HOME/ECMWFDATA$VERSION
# export PYTHONPATH=$ECMWFDATA/python
- export PATH=${PATH}:${HOME}/ECMWFDATA7.0/python
+ export PATH=${PATH}:${HOME}/ECMWFDATA7.1/python
;;
*cca*)
module switch PrgEnv-cray PrgEnv-intel
@@ -42,12 +42,12 @@ case $HOST in
module load emos
module load python
export SCRATCH=$TMPDIR
-# export ECMWFDATA=$HOME/ECMWFDATA7.0
+# export ECMWFDATA=$HOME/ECMWFDATA$VERSION
# export PYTHONPATH=$ECMWFDATA/python
- export PATH=${PATH}:${HOME}/ECMWFDATA7.0/python
+ export PATH=${PATH}:${HOME}/ECMWFDATA7.1/python
;;
# *)
-# export ECMWFDATA=$HOME/ECMWFDATA7.0
+# export ECMWFDATA=$HOME/ECMWFDATA$VERSION
# export PATH=/opt/anaconda/bin:$ECMWFDATA/python:${PATH}
# export PYTHONPATH=/opt/anaconda/lib/python2.7/site-packages/grib_api:$ECMWFDATA/python
# export SCRATCH=$ECMWFDATA/python
@@ -62,8 +62,8 @@ cd python$$
export CONTROL=CONTROL
cat >$CONTROL<<EOF
-GATEWAY srvx7.img.univie.ac.at
-DESTINATION leo@genericSftp
+GATEWAY srvx8.img.univie.ac.at
+DESTINATION philipa8@genericSftp
accuracy 16
addpar 186 187 188 235 139 39
basetime None
@@ -73,9 +73,9 @@ debug 0
dpdeta 1
dtime 3
ecfsdir ectmp:/${USER}/econdemand/
-ecstorage 1
+ecstorage 0
ectrans 0
-end_date 20131107
+end_date 20160809
eta 0
etadiff 0
expver 1
@@ -84,10 +84,10 @@ gauss 1
grib2flexpart 0
grid 5000
inputdir ../work
-left -175000
+left -15000
level 60
-levelist 1/to/60
-lower -90000
+levelist 55/to/60
+lower 30000
mailfail ${USER}
mailops ${USER}
makefile None
@@ -97,20 +97,20 @@ number OFF
omega 0
omegadiff 0
outputdir ../work
-prefix EN
+prefix EI
resol 63
-right 180000
+right 45000
smooth 0
-start_date 20131107
+start_date 20160809
step 00 01 02 03 04 05 00 07 08 09 10 11 00 01 02 03 04 05 00 07 08 09 10 11
stream OPER
time 00 00 00 00 00 00 06 00 00 00 00 00 12 12 12 12 12 12 18 12 12 12 12 12
type AN FC FC FC FC FC AN FC FC FC FC FC AN FC FC FC FC FC AN FC FC FC FC FC
-upper 90000
+upper 75000
EOF
-submit.py --controlfile=$CONTROL --inputdir=./work --outputdir=./work >prot
+submit.py --controlfile=$CONTROL --inputdir=./work --outputdir=./work 1> prot 2>&1
if [ $? -eq 0 ] ; then
l=0
diff --git a/python/job.temp b/python/job.temp
index d8291d0..c9cff9f 100644
--- a/python/job.temp
+++ b/python/job.temp
@@ -4,7 +4,7 @@
# start with ecaccess-job-submit -queueName ecgb NAME_OF_THIS_FILE on gateway server
# start with sbatch NAME_OF_THIS_FILE directly on machine
-#SBATCH --workdir=/scratch/ms/spatlh00/lh0
+#SBATCH --workdir=/scratch/ms/at/km4a
#SBATCH --qos=normal
#SBATCH --job-name=flex_ecmwf
#SBATCH --output=flex_ecmwf.%j.out
@@ -24,7 +24,7 @@
##PBS -l EC_memory_per_task=32000MB
set -x
-
+export VERSION=7.1
case $HOST in
*ecg*)
module load python
@@ -32,9 +32,9 @@ case $HOST in
module unload emos
module load grib_api/1.14.5
module load emos/437-r64
-# export ECMWFDATA=$HOME/ECMWFDATA7.0
+# export ECMWFDATA=$HOME/ECMWFDATA$VERSION
# export PYTHONPATH=$ECMWFDATA/python
- export PATH=${PATH}:${HOME}/ECMWFDATA7.0/python
+ export PATH=${PATH}:${HOME}/ECMWFDATA7.1/python
;;
*cca*)
module switch PrgEnv-cray PrgEnv-intel
@@ -42,12 +42,12 @@ case $HOST in
module load emos
module load python
export SCRATCH=$TMPDIR
-# export ECMWFDATA=$HOME/ECMWFDATA7.0
+# export ECMWFDATA=$HOME/ECMWFDATA$VERSION
# export PYTHONPATH=$ECMWFDATA/python
- export PATH=${PATH}:${HOME}/ECMWFDATA7.0/python
+ export PATH=${PATH}:${HOME}/ECMWFDATA7.1/python
;;
# *)
-# export ECMWFDATA=$HOME/ECMWFDATA7.0
+# export ECMWFDATA=$HOME/ECMWFDATA$VERSION
# export PATH=/opt/anaconda/bin:$ECMWFDATA/python:${PATH}
# export PYTHONPATH=/opt/anaconda/lib/python2.7/site-packages/grib_api:$ECMWFDATA/python
# export SCRATCH=$ECMWFDATA/python
@@ -62,14 +62,14 @@ cd python$$
export CONTROL=CONTROL
cat >$CONTROL<<EOF
-GATEWAY srvx7.img.univie.ac.at
-DESTINATION leo@genericSftp
+GATEWAY srvx8.img.univie.ac.at
+DESTINATION philipa8@genericSftp
EOF
cat >>$CONTROL<<EOF
EOF
-submit.py --controlfile=$CONTROL --inputdir=./work --outputdir=./work >prot
+submit.py --controlfile=$CONTROL --inputdir=./work --outputdir=./work 1> prot 2>&1
if [ $? -eq 0 ] ; then
l=0
diff --git a/python/job.temp.o b/python/job.temp.o
index cf59060..b3ff36c 100644
--- a/python/job.temp.o
+++ b/python/job.temp.o
@@ -24,7 +24,7 @@
##PBS -l EC_memory_per_task=32000MB
set -x
-
+export VERSION=7.1
case $HOST in
*ecg*)
module load python
@@ -32,7 +32,7 @@ case $HOST in
module unload emos
module load grib_api/1.14.5
module load emos/437-r64
-# export ECMWFDATA=$HOME/ECMWFDATA7.0
+# export ECMWFDATA=$HOME/ECMWFDATA$VERSION
# export PYTHONPATH=$ECMWFDATA/python
export PATH=${PATH}:
;;
@@ -42,12 +42,12 @@ case $HOST in
module load emos
module load python
export SCRATCH=$TMPDIR
-# export ECMWFDATA=$HOME/ECMWFDATA7.0
+# export ECMWFDATA=$HOME/ECMWFDATA$VERSION
# export PYTHONPATH=$ECMWFDATA/python
export PATH=${PATH}:
;;
# *)
-# export ECMWFDATA=$HOME/ECMWFDATA7.0
+# export ECMWFDATA=$HOME/ECMWFDATA$VERSION
# export PATH=/opt/anaconda/bin:$ECMWFDATA/python:${PATH}
# export PYTHONPATH=/opt/anaconda/lib/python2.7/site-packages/grib_api:$ECMWFDATA/python
# export SCRATCH=$ECMWFDATA/python
@@ -65,7 +65,7 @@ cat >>$CONTROL<<EOF
EOF
-submit.py --controlfile=$CONTROL --inputdir=./work --outputdir=./work >prot
+submit.py --controlfile=$CONTROL --inputdir=./work --outputdir=./work 1> prot 2>&1
if [ $? -eq 0 ] ; then
l=0
diff --git a/python/joboper.ksh b/python/joboper.ksh
index 3f1c284..94ab09f 100644
--- a/python/joboper.ksh
+++ b/python/joboper.ksh
@@ -4,7 +4,7 @@
# start with ecaccess-job-submit -queueName ecgb NAME_OF_THIS_FILE on gateway server
# start with sbatch NAME_OF_THIS_FILE directly on machine
-#SBATCH --workdir=/scratch/ms/spatlh00/lh0
+#SBATCH --workdir=/scratch/ms/at/km4a
#SBATCH --qos=normal
#SBATCH --job-name=flex_ecmwf
#SBATCH --output=flex_ecmwf.%j.out
@@ -24,7 +24,7 @@
##PBS -l EC_memory_per_task=32000MB
set -x
-
+export VERSION=7.1
case $HOST in
*ecg*)
module load python
@@ -32,9 +32,9 @@ case $HOST in
module unload emos
module load grib_api/1.14.5
module load emos/437-r64
-# export ECMWFDATA=$HOME/ECMWFDATA7.0
+# export ECMWFDATA=$HOME/ECMWFDATA$VERSION
# export PYTHONPATH=$ECMWFDATA/python
- export PATH=${PATH}:${HOME}/ECMWFDATA7.0/python
+ export PATH=${PATH}:${HOME}/ECMWFDATA7.1/python
;;
*cca*)
module switch PrgEnv-cray PrgEnv-intel
@@ -42,12 +42,12 @@ case $HOST in
module load emos
module load python
export SCRATCH=$TMPDIR
-# export ECMWFDATA=$HOME/ECMWFDATA7.0
+# export ECMWFDATA=$HOME/ECMWFDATA$VERSION
# export PYTHONPATH=$ECMWFDATA/python
- export PATH=${PATH}:${HOME}/ECMWFDATA7.0/python
+ export PATH=${PATH}:${HOME}/ECMWFDATA7.1/python
;;
# *)
-# export ECMWFDATA=$HOME/ECMWFDATA7.0
+# export ECMWFDATA=$HOME/ECMWFDATA$VERSION
# export PATH=/opt/anaconda/bin:$ECMWFDATA/python:${PATH}
# export PYTHONPATH=/opt/anaconda/lib/python2.7/site-packages/grib_api:$ECMWFDATA/python
# export SCRATCH=$ECMWFDATA/python
@@ -62,8 +62,8 @@ cd python$$
export CONTROL=CONTROL
cat >$CONTROL<<EOF
-GATEWAY srvx7.img.univie.ac.at
-DESTINATION leo@genericSftp
+GATEWAY srvx8.img.univie.ac.at
+DESTINATION philipa8@genericSftp
accuracy 16
addpar 186 187 188 235 139 39
basetime None
@@ -73,7 +73,7 @@ debug 0
dpdeta 1
dtime 3
ecfsdir ectmp:/${USER}/econdemand/
-ecstorage 1
+ecstorage 0
ectrans 0
start_date ${MSJ_YEAR}${MSJ_MONTH}${MSJ_DAY}
eta 0
@@ -84,10 +84,10 @@ gauss 1
grib2flexpart 0
grid 5000
inputdir ../work
-left -175000
+left -15000
level 60
-levelist 1/to/60
-lower -90000
+levelist 55/to/60
+lower 30000
mailfail ${USER}
mailops ${USER}
makefile None
@@ -97,20 +97,20 @@ number OFF
omega 0
omegadiff 0
outputdir ../work
-prefix EN
+prefix EI
resol 63
-right 180000
+right 45000
smooth 0
start_date ${MSJ_YEAR}${MSJ_MONTH}${MSJ_DAY}
step 00 01 02 03 04 05 00 07 08 09 10 11 00 01 02 03 04 05 00 07 08 09 10 11
stream OPER
time 00 00 00 00 00 00 06 00 00 00 00 00 12 12 12 12 12 12 18 12 12 12 12 12
type AN FC FC FC FC FC AN FC FC FC FC FC AN FC FC FC FC FC AN FC FC FC FC FC
-upper 90000
+upper 75000
EOF
-submit.py --controlfile=$CONTROL --inputdir=./work --outputdir=./work >prot
+submit.py --controlfile=$CONTROL --inputdir=./work --outputdir=./work 1> prot 2>&1
if [ $? -eq 0 ] ; then
l=0
diff --git a/python/plot_retrieved.py b/python/plot_retrieved.py
old mode 100644
new mode 100755
index b95df2d..7d4c6ad
--- a/python/plot_retrieved.py
+++ b/python/plot_retrieved.py
@@ -29,7 +29,8 @@ import matplotlib.cm as cmx
import matplotlib.colors as colors
from argparse import ArgumentParser,ArgumentDefaultsHelpFormatter
-from FlexpartTools import interpret_args_and_control,silentremove,product,Control
+from Tools import interpret_args_and_control, silentremove, product
+from Control import Control
from GribTools import GribTools
from gribapi import *
from rasotools.utils import stats
diff --git a/python/prepareFLEXPART.py b/python/prepareFLEXPART.py
old mode 100644
new mode 100755
index 611b56b..9c5ccbe
--- a/python/prepareFLEXPART.py
+++ b/python/prepareFLEXPART.py
@@ -29,9 +29,11 @@
- applied PEP8 style guide
- added documentation
- minor changes in programming style for consistence
+ - BUG: removed call of cleanup-Function after call of prepareFlexpart
+ since it is already called in prepareFlexpart at the end!
@License:
- (C) Copyright 2014.
+ (C) Copyright 2014-2018.
This software is licensed under the terms of the Apache Licence Version 2.0
which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
@@ -39,7 +41,6 @@
@Requirements:
- A standard python 2.6 or 2.7 installation
- dateutils
- - matplotlib (optional, for debugging)
- ECMWF specific packages, all available from https://software.ecmwf.int/
ECMWF WebMARS, gribAPI with python enabled, emoslib and
ecaccess web toolkit
@@ -51,33 +52,38 @@
Prepare input 3D-wind fields in hybrid coordinates +
surface fields for FLEXPART runs
"""
-
+# ------------------------------------------------------------------------------
+# MODULES
+# ------------------------------------------------------------------------------
import calendar
import shutil
import datetime
import time
-import os,inspect,sys
+import os
+import inspect
+import sys
import socket
from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter
-# add path to submit.py to pythonpath so that python finds its buddies
-localpythonpath=os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
-if localpythonpath not in sys.path:
- sys.path.append(localpythonpath)
-from UIOTools import UIOFiles
-#from string import strip
-from GribTools import GribTools
-from FlexpartTools import ECFlexpart, Control, interpret_args_and_control, cleanup
-
-hostname=socket.gethostname()
-ecapi= 'ecmwf' not in hostname
+import UIOFiles
+import Control
+import Tools
+import ECFlexpart
+hostname = socket.gethostname()
+ecapi = 'ecmwf' not in hostname
try:
if ecapi:
import ecmwfapi
except ImportError:
ecapi = False
-
+# add path to submit.py to pythonpath so that python finds its buddies
+localpythonpath=os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
+if localpythonpath not in sys.path:
+ sys.path.append(localpythonpath)
+# ------------------------------------------------------------------------------
+# FUNCTION
+# ------------------------------------------------------------------------------
def prepareFLEXPART(args, c):
'''
@Description:
@@ -120,14 +126,14 @@ def prepareFLEXPART(args, c):
month=int(c.end_date[4:6]),
day=int(c.end_date[6:]))
- # to deaccumulated the fluxes correctly
+ # to deaccumulate the fluxes correctly
# one day ahead of the start date and
# one day after the end date is needed
startm1 = start - datetime.timedelta(days=1)
endp1 = end + datetime.timedelta(days=1)
# get all files with flux data to be deaccumulated
- inputfiles = UIOFiles(['.grib', '.grb', '.grib1',
+ inputfiles = UIOFiles.UIOFiles(['.grib', '.grb', '.grib1',
'.grib2', '.grb1', '.grb2'])
inputfiles.listFiles(c.inputdir, '*OG_acc_SL*.' + c.ppid + '.*')
@@ -137,7 +143,7 @@ def prepareFLEXPART(args, c):
os.makedirs(c.outputdir)
# deaccumulate the flux data
- flexpart = ECFlexpart(c, fluxes=True)
+ flexpart = ECFlexpart.ECFlexpart(c, fluxes=True)
flexpart.write_namelist(c, 'fort.4')
flexpart.deacc_fluxes(inputfiles, c)
@@ -145,7 +151,7 @@ def prepareFLEXPART(args, c):
"/to/" + end.strftime("%Y%m%d"))
# get a list of all files from the root inputdir
- inputfiles = UIOFiles(['.grib', '.grb', '.grib1',
+ inputfiles = UIOFiles.UIOFiles(['.grib', '.grb', '.grib1',
'.grib2', '.grb1', '.grb2'])
inputfiles.listFiles(c.inputdir, '????__??.*' + c.ppid + '.*')
@@ -156,20 +162,19 @@ def prepareFLEXPART(args, c):
if c.basetime == '00':
start = startm1
- flexpart = ECFlexpart(c, fluxes=False)
+ flexpart = ECFlexpart.ECFlexpart(c, fluxes=False)
flexpart.create(inputfiles, c)
flexpart.process_output(c)
# check if in debugging mode, then store all files
+ # otherwise delete temporary files
if int(c.debug) != 0:
print('Temporary files left intact')
else:
- cleanup(c)
+ Tools.cleanup(c)
return
if __name__ == "__main__":
- args, c = interpret_args_and_control()
+ args, c = Tools.interpret_args_and_control()
prepareFLEXPART(args, c)
- cleanup(c)
-
diff --git a/python/profiling.py b/python/profiling.py
new file mode 100644
index 0000000..c03d35d
--- /dev/null
+++ b/python/profiling.py
@@ -0,0 +1,69 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#************************************************************************
+# TODO AP
+# - add description of file
+# - check of license of book content
+#************************************************************************
+"""
+@Author: Anne Philipp (University of Vienna)
+
+@Date: March 2018
+
+@License:
+ (C) Copyright 2018.
+
+ This software is licensed under the terms of the Apache Licence Version 2.0
+ which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+
+@Requirements:
+ A standard python 2.6 or 2.7 installation
+
+@Description:
+ ...
+"""
+# ------------------------------------------------------------------------------
+# MODULES
+# ------------------------------------------------------------------------------
+from functools import wraps
+import time
+
+# ------------------------------------------------------------------------------
+# FUNCTION
+# ------------------------------------------------------------------------------
+def timefn(fn):
+ '''
+ @Description:
+ Decorator function. It takes the inner function as an argument.
+ '''
+ @wraps(fn)
+ def measure_time(*args, **kwargs):
+ '''
+ @Descripton:
+ Passes the arguments through fn for execution. Around the
+ execution of fn the time is captured to execute the fn function
+ and prints the result along with the function name.
+
+ This is taken from the book "High Performance Python" from
+ Micha Gorelick and Ian Ozsvald, O'Reilly publisher, 2014,
+ ISBN: 978-1-449-36159-4
+
+ @Input:
+ *args: undefined
+ A variable number of positional arguments.
+
+ **kwargs: undefined
+ A variable number of key/value arguments.
+
+ @Return:
+ <nothing>
+ '''
+
+ t1 = time.time()
+ result = fn(*args, **kwargs)
+ t2 = time.time()
+ print("@timefn:" + fn.func_name + " took " + str(t2 - t1) + " seconds")
+
+ return result
+
+ return measure_time
diff --git a/python/submit.py b/python/submit.py
old mode 100644
new mode 100755
index e3dcb59..3c4030c
--- a/python/submit.py
+++ b/python/submit.py
@@ -2,10 +2,13 @@
# -*- coding: utf-8 -*-
#************************************************************************
# TODO AP
-#AP
+#
# - Change History ist nicht angepasst ans File! Original geben lassen
# - dead code ? what to do?
# - seperate operational and reanlysis for clarification
+# - add correct file description
+# - divide in two submits , ondemand und operational
+# -
#************************************************************************
"""
@Author: Anne Fouilloux (University of Oslo)
@@ -14,59 +17,42 @@
@ChangeHistory:
November 2015 - Leopold Haimberger (University of Vienna):
- - using the WebAPI also for general MARS retrievals
- job submission on ecgate and cca
- job templates suitable for twice daily operational dissemination
- - dividing retrievals of longer periods into digestable chunks
- - retrieve also longer term forecasts, not only analyses and
- short term forecast data
- - conversion into GRIB2
- - conversion into .fp format for faster execution of FLEXPART
February 2018 - Anne Philipp (University of Vienna):
- applied PEP8 style guide
- added documentation
- - minor changes in programming style for consistence
+ - minor changes in programming style (for consistence)
@License:
- (C) Copyright 2014 UIO.
+ (C) Copyright 2014-2018.
This software is licensed under the terms of the Apache Licence Version 2.0
which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
@Requirements:
- A standard python 2.6 or 2.7 installation
- - dateutils
- - matplotlib (optional, for debugging)
- - ECMWF specific packages, all available from https://software.ecmwf.int/
- ECMWF WebMARS, gribAPI with python enabled, emoslib and
- ecaccess web toolkit
@Description:
Further documentation may be obtained from www.flexpart.eu.
- Functionality provided:
- Prepare input 3D-wind fields in hybrid coordinates +
- surface fields for FLEXPART runs
+
"""
# ------------------------------------------------------------------------------
# MODULES
# ------------------------------------------------------------------------------
-import calendar
-import shutil
-import datetime
-import time
-import os, sys, glob
+import os
+import sys
+import glob
import subprocess
import inspect
-# add path to submit.py to pythonpath so that python finds its buddies
+# add the pythondir path so that python finds its buddies (from flex_extract)
localpythonpath = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
sys.path.append(localpythonpath)
-from UIOTools import UIOFiles
-from string import strip
-from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter
-from GribTools import GribTools
-from FlexpartTools import ECFlexpart, Control, interpret_args_and_control, normalexit, myerror
+
+# software specific classes and modules from flex_extract
+from Tools import interpret_args_and_control, normalexit
from getMARSdata import getMARSdata
from prepareFLEXPART import prepareFLEXPART
# ------------------------------------------------------------------------------
@@ -75,9 +61,10 @@ from prepareFLEXPART import prepareFLEXPART
def main():
'''
@Description:
- Get the arguments from script call and initialize an object from
- Control class. Decides from the argument "queue" if the local version
- is done "queue=None" or the gateway version "queue=ecgate".
+ Get the arguments from script call and from Control file.
+ Decides from the argument "queue" if the local version
+ is done "queue=None" or the gateway version with "queue=ecgate"
+ or "queue=cca".
@Input:
<nothing>
@@ -87,6 +74,7 @@ def main():
'''
calledfromdir = os.getcwd()
args, c = interpret_args_and_control()
+ # on local side
if args.queue is None:
if c.inputdir[0] != '/':
c.inputdir = os.path.join(calledfromdir, c.inputdir)
@@ -95,13 +83,13 @@ def main():
getMARSdata(args, c)
prepareFLEXPART(args, c)
normalexit(c)
+ # on ECMWF server
else:
submit(args.job_template, c, args.queue)
return
def submit(jtemplate, c, queue):
- #AP divide in two submits , ondemand und operational
'''
@Description:
Prepares the job script and submit it to the specified queue.
@@ -137,15 +125,13 @@ def submit(jtemplate, c, queue):
with open(jtemplate) as f:
lftext = f.read().split('\n')
insert_point = lftext.index('EOF')
-#AP es gibt mehrere EOFs überprüfen!
# put all parameters of control instance into a list
- clist = c.tolist() # reanalysis (EI)
+ clist = c.tolist()
colist = [] # operational
mt = 0
#AP wieso 2 for loops?
-#AP dieser part ist für die CTBTO Operational retrieves bis zum aktuellen Tag.
for elem in clist:
if 'maxstep' in elem:
mt = int(elem.split(' ')[1])
@@ -161,9 +147,6 @@ def submit(jtemplate, c, queue):
if 'time' in elem and mt > 24:
elem = 'time ' + '${MSJ_BASETIME} {MSJ_BASETIME}'
colist.append(elem)
-#AP end
-
-#AP whats the difference between clist and colist ?! What is MSJ?
lftextondemand = lftext[:insert_point] + clist + lftext[insert_point + 2:]
lftextoper = lftext[:insert_point] + colist + lftext[insert_point + 2:]
@@ -171,17 +154,13 @@ def submit(jtemplate, c, queue):
with open('job.ksh', 'w') as h:
h.write('\n'.join(lftextondemand))
-#AP this is not used ?! what is it for?
-#maybe a differentiation is needed
- h = open('joboper.ksh', 'w')
- h.write('\n'.join(lftextoper))
- h.close()
-#AP end
+ with open('joboper.ksh', 'w') as h:
+ h.write('\n'.join(lftextoper))
# submit job script to queue
try:
p = subprocess.check_call(['ecaccess-job-submit', '-queueName',
- queue,' job.ksh'])
+ queue, 'job.ksh'])
except:
print('ecaccess-job-submit failed, probably eccert has expired')
exit(1)
diff --git a/python/testsuite.py b/python/testsuite.py
old mode 100644
new mode 100755
diff --git a/pythontest/TestTools.py b/pythontest/TestTools.py
new file mode 100644
index 0000000..5fe4feb
--- /dev/null
+++ b/pythontest/TestTools.py
@@ -0,0 +1,37 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import unittest
+import sys
+sys.path.append('../python')
+from Tools import init128, toparamId
+
+
+class TestTools(unittest.TestCase):
+ '''
+ '''
+
+ def setUp(self):
+ pass
+
+ def test_init128(self):
+ '''
+ '''
+ table128 = init128('../grib_templates/ecmwf_grib1_table_128')
+ expected = {'078': 'TCLW', '130': 'T', '034': 'SST'}
+ # check a sample of parameters which must have been read in
+ result = all((k in table128 and table128[k]==v) for k,v in expected.iteritems())
+ self.assertEqual(result, True)
+
+
+ def test_toparamId(self):
+ '''
+ '''
+ table128 = init128('../grib_templates/ecmwf_grib1_table_128')
+ pars = toparamId("T/SP/LSP/SSHF", table128)
+ for par in pars:
+ self.assertIn(par, [130, 134, 142, 146])
+
+
+if __name__ == "__main__":
+ unittest.main()
\ No newline at end of file
diff --git a/pythontest/TestUIOFiles.py b/pythontest/TestUIOFiles.py
new file mode 100644
index 0000000..3d929cb
--- /dev/null
+++ b/pythontest/TestUIOFiles.py
@@ -0,0 +1,64 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import unittest
+import os
+import sys
+sys.path.append('../python')
+import UIOFiles
+
+
+class TestUIOFiles(unittest.TestCase):
+ '''
+ Test class to test the UIOFiles methods.
+ '''
+
+ def setUp(self):
+ '''
+ @Description:
+ Prepare test case. Initialize comparing filelist and
+ the test path.
+
+ @Input:
+ self: instance of TestUIOFiles
+ Class to test the UIOFiles methods.
+
+ @Return:
+ <nothing>
+ '''
+ self.testpath = os.path.join(os.path.dirname(__file__), 'TestDir')
+ self.expected = ['FCGG__SL.20160410.40429.16424.grb',
+ 'FCOG__ML.20160410.40429.16424.grb',
+ 'FCSH__ML.20160410.40429.16424.grb',
+ 'OG_OROLSM__SL.20160410.40429.16424.grb',
+ 'FCOG_acc_SL.20160409.40429.16424.grb',
+ 'FCOG__SL.20160410.40429.16424.grb',
+ 'FCSH__SL.20160410.40429.16424.grb']
+
+ return
+
+ def test_listFiles(self):
+ '''
+ @Description:
+ Test the listFiles method from class UIOFiles.
+
+ @Input:
+ self: instance of TestClass
+ Class to test the UIOFiles methods.
+
+ @Return:
+ <nothing>
+ '''
+
+ # Initialise and collect filenames
+ files = UIOFiles.UIOFiles(['.grb'])
+ files.listFiles(self.testpath, '*')
+ # get the basename to just check for equality of filenames
+ filelist = [os.path.basename(f) for f in files.files]
+ # comparison of expected filenames against the collected ones
+ self.assertItemsEqual(self.expected, filelist)
+
+ return
+
+if __name__ == "__main__":
+ unittest.main()
\ No newline at end of file
diff --git a/pythontest/__init__.py b/pythontest/__init__.py
new file mode 100644
index 0000000..749a109
--- /dev/null
+++ b/pythontest/__init__.py
@@ -0,0 +1,13 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""
+@Author: Anne Philipp (University of Vienna)
+
+@Date: March 2018
+
+@License:
+ (C) Copyright 2014 UIO.
+
+ This software is licensed under the terms of the Apache Licence Version 2.0
+ which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+"""
\ No newline at end of file
diff --git a/src/Makefile.local.ifort b/src/Makefile.local.ifort
index 206863f..ef6a78b 100644
--- a/src/Makefile.local.ifort
+++ b/src/Makefile.local.ifort
@@ -11,8 +11,8 @@
.s .s~ .sh .sh~ .h .h~ .C .C~ .a
-GRIB_API_INCLUDE_DIR=/usr/local/ifort/grib_api-1.14.3//include
-GRIB_API_LIB=-openmp -L/usr/local/ifort/grib_api-1.14.3/lib -Bstatic -lgrib_api_f77 -lgrib_api_f90 -lgrib_api -Bdynamic -lm -ljasper #-lopenjpeg
+GRIB_API_INCLUDE_DIR=/usr/local/ifort/grib1.12.3//include
+GRIB_API_LIB=-openmp -L/usr/local/ifort/grib1.12.3//lib -Bstatic -lgrib_api_f77 -lgrib_api_f90 -lgrib_api -Bdynamic -lm -ljasper #-lopenjpeg
OPT = -g
DEBUG = -g
--
GitLab