From 4b644c9d3c6aef257f13244dd381b73f1a0e3ae8 Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Thu, 11 Jan 2024 10:31:54 -0800 Subject: [PATCH 1/2] add new data sets for WACCM --- data/cross_sections/BRO_JPL06.nc | Bin 0 -> 11768 bytes data/cross_sections/CL2O2_JPL10.nc | Bin 0 -> 17686 bytes data/cross_sections/CLO_JPL06.nc | Bin 0 -> 9189 bytes data/cross_sections/HNO3_JPL06.nc | Bin 0 -> 3436 bytes examples/ts1_tsmlt.json | 24 ++- test/data/xsqy.doug.config.json | 121 +++++++++++++ test/unit/tuv_doug/JCALC/CMakeLists.txt | 4 + test/unit/tuv_doug/JCALC/XSQY_BRO.f | 116 +++++++++++++ test/unit/tuv_doug/JCALC/XSQY_CL2O2.f | 117 +++++++++++++ test/unit/tuv_doug/JCALC/XSQY_CLO.f | 115 +++++++++++++ test/unit/tuv_doug/JCALC/XSQY_HNO3.f | 140 +++++++++++++++ test/unit/tuv_doug/data_sets.F90 | 25 ++- test/unit/tuv_doug/driver.F90 | 12 ++ tool/data_conversion/photo.config.json | 44 +++++ tool/data_conversion/text_to_netcdf.py | 61 +++++++ tool/data_conversion/xsqy_subs.py | 218 ++++++++++++++++++++++++ 16 files changed, 988 insertions(+), 9 deletions(-) create mode 100644 data/cross_sections/BRO_JPL06.nc create mode 100644 data/cross_sections/CL2O2_JPL10.nc create mode 100644 data/cross_sections/CLO_JPL06.nc create mode 100644 data/cross_sections/HNO3_JPL06.nc create mode 100644 test/unit/tuv_doug/JCALC/XSQY_BRO.f create mode 100644 test/unit/tuv_doug/JCALC/XSQY_CL2O2.f create mode 100644 test/unit/tuv_doug/JCALC/XSQY_CLO.f create mode 100644 test/unit/tuv_doug/JCALC/XSQY_HNO3.f create mode 100644 tool/data_conversion/photo.config.json create mode 100644 tool/data_conversion/text_to_netcdf.py create mode 100644 tool/data_conversion/xsqy_subs.py diff --git a/data/cross_sections/BRO_JPL06.nc b/data/cross_sections/BRO_JPL06.nc new file mode 100644 index 0000000000000000000000000000000000000000..5ad156c8a4f773e3b51b541d4c920e4dace4e252 GIT binary patch literal 11768 zcmeI2d0dp`zQA7?)^Su4Qxun1DMcN|VOSjChYicHFXMoSRxkrIFe9_j%sU8`g=J=n zMwa^$WGYtnp>wWVTA`I-rI(h|#U;ZfB{e1V6vs->^UV8u#bf8*dpw_W?&to&gUpxr zS%2H}d*|cxypWchGEh8H?CIyn2;~-*-Y~vKjtz&lEt;8|maO%Z{DINoNcQ zpwBBr^pj!Ic<@a^V+8F>N*)j9gU5%jb&Eex>ESVyR(g7PcpSo3fuf;w{`w5y2mVf^UaNm&Vbc}yUwqN-3?lhd=+ zIeF>1ImLNN30Z1JGKVT8V!ZhQISJWncYhW3&o`N!EcvHrZALCBIfXTttwx*EY`3u` zE{?U?xxS_36PdZbbs*DaQP`pA)WF)i4N6KPka;eEs78M^Z z3s=Y#(++PH@or)6aEm)?$asPOm|GZT8oFf#8B2@c7O#JuTlmmBj-hAb*SR-iy0sdF zBkw@){QPIAdHE%BaCc=0J$PgIjHb^XUV>@Fm3oKXYUGHe7g3Xwvc!M(Z*DjJ8zYEi z!ydJ5Tq>PEM`;6Bmq$%GyPX%OGjuHtRG25tLK zxzBCo0d(Cv%5U2q()IaH;>$PPm6?3T(rOPc^b;fDj8?L3^qk9KbSCmij`{KF1xNq= z-jg$Rb-sK)A9&}^6B407^d9>V+3K{Yn8*Q4T(1 zOGs#m$Pd0=TsAYur(zmaG~PbU>9-Tb19Y?})tUZ5iV$jdKi^yE@jjNjXJ?o-=p-TV z{)ZO;d6(T|;8CS$5hw`Gi*N*^XD;7~VelCro6kNtuP@XpFpu02I`eb+xMS4Ai}we^ z1eH4MPG_;xSjv&tX>ni5ufUjY%=13vcE%(38hC;d0^I^&;87#q`Ao$*abR1=oq;H1 zVnI}ig%8BRjGTtPR-K~Os&kUm#qQ9LCd<%3Gw8X0IuHe&1(RKjNqzqBNfd-JAbFAo zHrDv1PJ3A}r22oef`@p0wS3@V`y-SH&<3zB=4m~GCh_#~0@`^cC7QkCH1Om>I+ z+Y~XQJJ#}N7O=_0aBN|j-C;F4*l?Mgj-vbCj~`X)Fp^wpa1WzzRVtxEu80hihlhpB z+3=`Xc~oqaQWC+N1otdw<}60{m>?XJMB|Occw_tGO){CL&D&o}%gi{A`joC!m{**r z%`MK*WXYptTs7C*ofSf=`^*1N3#2CHX1Rlu`i~SLV1;+O5Sv1v3-!V`P26Gbw-_BS z)C&=(j+>1^O+}Z-IUb^jo0c0{s@~w?MxI`Yq6Jf&U{5 zbUmhWuhfm~MKTV)-VYgo3_`NVP^28GM5>T!NDZW>US1|eBwC{m78B2~yVqy|}t)FDmCN@O*%7P%N% zk6eRnKsF(pk$aFW$Tnm_T=Udyx!>`Xd97K}Z%Eij*UjNEI>-sX-PZbx0Gk z5?PI`MJ`6xBiA4skWI*D_B!RyO7<;UL@l}{gDC4AS8EUC3@^FOsQ3{gDC4AS8EUC3@^FOr#u z`Xd97K}Z%Eij*UjNEI>-sX-PZbx0Gk5?PI`MJ`6xBiA4skWI*D_B!R zyO7<;UL;eE`Xd9XRrC;pKjnHGk67KOXAHwg93pZ6AR-4X?&~ck1AlZeoBR(k{moAs*(_JZRb6HEsX@!0Yz73Yv^ z4zd0v4!hNCGqP5@&0cLb7{&ZSN$6~&-Z3fGpPW2|Ddh6#Q2H@RoZolA9_z1&QD%x~ z-9EHs9eP_inNgvPmPLg{D5GR@n>CaOHA$RkF0q*H;+YZP`;P0@j4cr`r+6!8IyepLmd_j!vN#+TzO=x8 zb$T4Uu)VD=G$9gp`ObwQGp51W^*&!6zor7+lS7iLmnq?wC8K&DZcxIaYl)I@V**$$ zm8CD<0th;?_RCM-mz*BtVbPC&wbr$AI&NlZsbAiHF2n+)?es zOlaIbYWBh_NighmP{p6)rbBtI*DFs)0@Sr^Y#Y5z1NPp)o;PdL;btAP?`%f`MCU&8 zOY1R!;)I^<2g5TU>)f!W>WVyAbjHh2VCd_YA5{3q`5RPV+x(664>Ap`O22>YsL8XS zYUgVgJc_cR$YaN&(F2pAeb1QxT7F+5e6oZqeqEgo55D(v&zVt=fRC<6eCg?#keRkk z`|)#`aCzawPrYnNhStc>_LR*^g7W8XJyF!Eg-7>#yj>KngKZ~j8a*z|f)xwCUGDQv zF06|STRHG%GJJk?LAdHf0<6B_DvyuX!q%f#mVR4N0@oKzZE=+qfz7XTPvUp^Q1A12 z(?{}5h;6OydG~N46l>+$jfzCb2})|2@Qw}wUXu?jXfuF%=Z%~UQ4tJ0ditooG9MOM zKAL*ztt_ZI)cNAcA~ig|**8$_TtgAVO97YgTfY|1#M1Loz`IXB}A;H$vP zQ#ZHe!nZ}9^EMWz!vtI50g*fv9#>qh|8!U;)Gn#MQ6E1SM)dx$^W8!QpsM@eg57&FA^PHM`?QcW7(8fV@Re6Gz+AOu)qq=hF!@1E$U{wKP^7AD zdupW>vh8UzXJ|@c%!3o+U;4ZN)~mFs=MH7Tg1BRC?`YKUesj&b=^K;bG3D&($}PEY z$)i5fM^OSj!TYq0fmSFOJ>QsbG{Lj;Dw5w!)4`NWjdPE5<-?ZazAf)%WP%!Mrqu6C zgzTHrFZb)G!*jj|H_!Pg6&{%Xhph#@MtJ3ml4q@7Tj9R>&LcCT%iyQ+AGO{ZS_EsG zf4LAjVkTUfw7H=9pN!nmtCXgYc68TGUI5L{F;^n>e}uzpzOOn*xXgeS?qS}`~ob|o(P zy7etJgnTxmBkH>h=s9`eLgCUpSgarYjb^S6HZ(SRp3#(HT$R8N+>5(^y;TU~gMakD z|BGCxKe>9>_kT);#z)^DcKxMzketz#E%k^4&yXuSc3g@FkF!nsajmIvKA~k{tBOU}NV+*$OY&!U?AF`$4L=q_W+Wc-F zQo-2CnODo#O@~X$-{M4m)8YA4-|q7xr-N>p=CXFV3UYIsTSK1Dg7j;Z%RY@Lf&cV< z-_Lfl0RF6~U%sj)9ft1I?po0_9p{(FIRppNrS@e?Dfj?$-rek zv3=bs6~yjJ%;;>3y*+dB%An6z$8RoREpU>p5;c=He?e6V!?oY(=;=x+sd1za3( z@}sFR@jjQIGAyEZ}IM?DIj_E_`xFY$&meaFqio2c-Wx4 z9B-X45}HbS4{v&67~Fc}@EA$bPzZc|XqErfA@K2r!5db8Hw0wP72Mk+0-!l&V7ikG!8lw_f1U-bO-zp*ZpDR3d>*4$Gw^QT*2Cd*W$cuR_%vz4Ez;9arvwO literal 0 HcmV?d00001 diff --git a/data/cross_sections/CL2O2_JPL10.nc b/data/cross_sections/CL2O2_JPL10.nc new file mode 100644 index 0000000000000000000000000000000000000000..1656cee42a64cb66c126ce6ba9fa650b2f0cd93d GIT binary patch literal 17686 zcmeI(d0Z3Mzc=t<5djeq7jPeODOd;rMA36BVG~5b2#C0n0D(k-q)9;BvDJ!-iYr*F zq7_Ap5?m@>tXTIAR;@_ZGEQ9A_`4Z~Wh}X`{1cs=_@5V8 z@n0Dx!i4D~^M`)?wcHVFg4sFVWvQ|IjRPz!n)4?uEi5eTaFtoB=KS-UcI@z>Nf%qj zWeRt71LIh%?D?jxn1lAF$F4T1lxagk`?4vD3|Xc| zq0D5Hvb1cbQv30vY^hoyO-h#uSkvL}>%%y39kS$}?}(5`%e2WUsog^r8H$vD3Od&& zcZ_~~ByeCzM1Kb^HTbSMn!1)49U@jMRWh|!A=4Q9=*C^lk11QG<~!8K*OTib_w+v9 z1-%49PkBnJe;*%#kFUVD=fSmBrfFgRG%d!g;qEv4Z<`i|>4|Asz}?Fq=ESr9&zlxg z>YYdG$xZG2y3`Ei2JzQhG0nMZ#8u}+RHG23Qk|m6lxn%c z@((Wpu6CPOfLZkPDS9j-5a)&v>4qR$kaewf7$;wn{O{JL6Dn8$mGe-Ho z>+-Wgs*!0wHvQq9!PxRWGNmBF?{1&V~-rTg2cZ^{;6fYRsLUQ#)&s` z`h7^X&ivh${FwP~xWk-kfBpIIrJiYyBLt>reuD`ySCIe93xIo<{msPV2gM&|%^f#& z;V80awCN&leCBVAHXAj1$j4M0gF0eqeuhaKV~$!_qa-k{$!euWlcVId z200z~rpI+-;*I#biDw>!Sz`x~czmY+`Dea67j~_y`Ew$gE3vt&4(X;uY{Cfn)Q<`4 z8x|859TJvkO#L2QGkl_@NVWg&L^O}t7*6iX-Cp4{{(+L-VA@%ox3i@1YI9;<@^7?zwf6B2)7SKRF0Z>!_jWm#UyFtz4-#{x-#h(U@y- zgX7syE^ur@s#2XHQ?ot-A%7RY?rlwXC97pztxPfA#=lkh9{LJ>yL$=LMVb_D1j=N3kzT|)WIsKhqcfEnbBw;vQP*`Pz)td1#@8mEQUH* z1@*8N8Xz+U?L!s{p$Lkh1gc;zEP%yO2dkhS)R~N3Kqdq2Llz352#TQus$ecGfW=S;tDqj%LIY$n(LQ9M5Q?A}N}vkn!U9+fb+8KR zVJ$R3Mv3+z3x!Yw#ZUrOFc%iUVyJ^vP!DUN0WvC6`?k;tx=+A{YV1Fab)S z9I9Xr%!PTd02aYwSPFHp5>~+jP!DThEv$nE_yRItV*Jnvx=+A{YV1Fab)S z9I9Xr%!PTd02aYwSPFHp5>~+jP!DThEv$nE_yRI&j2}8dSI9zlD1-q}1S6mrCO`?4 zLlw+{xiAkFz#>=-OQ8-{!YX(G>R}D6g>}#XUqD8K@k1x*3R&n5g)jh$U<4Gy1So-W zsDe2#7v{kNSOkk2KnavX z70iLTFb@{MB3KMdp$=BUDtG|uVGXQ>bcB~T7k zFbC$sJXipWU@y--~p(IHLw=eK?8gNneiAubb_vsh3-%Y1E2_sp#-X6E-Zk> zPzS4^9@atwWG0yU=LA_Ogd!-05~zZ?umBcA9jtd6D)94TS#C2jw=3DM6Pv0~YqTuCHHcNqlDUn4o@{1TMiRHtMVb2N(;C)3R;$iR z=C@7uwhKugI4~sLKAPJE+uJThkt|PFsuh_sZacIxN0B1)WIvqAP-YrW8Fxyn|MUU3 z2h;e0y=l{p`%sxw?a|v#ELBM}r5bLhl0uo~?kg1b^01HlxMQ=o9e1j~{g4m409my( zGmYy6ciOAFFSql>K3I{Iu280_r7HQj-uCvPg8l;bAD-&T#tTAB=eVtS3U0@{GRu<< z6$}xueWh~QI8RoZnZia2f(2}FmU@g_nVur!dhW*U1d=9clxh{fAB@$gxQ$`ztPDO* zTDgqnLSh;t>mHKM?Z)w7lVsYlGFc|ebIx zA6db5C|#PNGB(JaR47wqUQ)GMI*z+Q%k@~PCfr18l<8TSH+Dj*9>+#;H*@>l)7@BA78iewoPSih zlIz-7g;s8ikd)i2o}4aY+WwlpU7`}`gk%?qev7AeFUnMRQvzliOJ2nm%l%knC`M1Etix^W0*bu7o_Q5 z(?@Gu3H@1Rl2*zUixf6#9Q!dQ{F=g2lI3v4f&Y&cT4}ne1{ssw{@B;<)xUVKeb`m&5?>A`=)-Qeye6bO9-d|2*^|Ka)U^&da%xf+m^B4yn(Gd%3Q(ry)tQ6n2Kl_EN!EktZu#{~J1WHV=MEJ^`^Spu$qo&= zu4BY>@y3!6Md&~}txXV_HY1vDYjfu$xfVrlwV9v(gho-9w!IqsP&j}_w7t4K>5u+& zT-%vFT1WM#_c!%XW$uim?aFW3c5NL=3xd849-rTjTAhlXSp6=7{`8{BkezuE)N}Tp zf@ZJ7>DBV-W2~o#)524I`#y6Fr*mgF|NSdnUphVbS|6vFzI4*d@}HJJ38Oy-&pNZe zEQ}gvr#yKq3!~|2{VOVagwflL!)S@FOZBx!p)|RDE3dvap|oSZL;m`zkDAY$ zuQT6oe%{MdOU9<`38NP6557L*)|c9+y`IudG&j7|<<(PAT{A zz4%zfhwpiHK?VoHrfnYsK_gNOg}d#}A@6^UEI2_&kPg&9{!%>tpEN_KH#4=ETx% zY2hu8d=W?Yz5Jnl#p^g~>hECmGoI4mTQzU`4WXIsj{E4G6VSf|S}}`V7*RBo{@!-_ zGOxkI=&4P2zIS#WPH&2$%j&C#)6+9OQygj+6lOiIa|6TSM`OA$W_jhoZerewQ>28JG5Htj13?e;puCnZWk zr_ZTgF~V8$q5r-ooMYaxiL@kT|LmBBqo`~BguQDUjH0t=3g5pPJCa^HF{&>2@(3!H zHFG;H9znm~6x3+=_Tltdz4B>oyWuoEXvL+}X~XE!Sx5HWzdMxPJNc?yJ!mLxdunw= zx3UBpbt>}eRhtC5vfa2tX&FQ4LM=;Cr@K&S;&dwBh@ifDqmO@{bmqtIVWyL z;5UQlk%|^cLZ?A=>X|dN}-6ryTAYIC~(q6%CuTM;U2sDV_f}d*B35xigF?`u732;EN$YNuzs{pJCCI=wntEN-lv{gv9ZtA zaC$B^_V)a)AJ?mS-I&*@d3~AJ;Thf0)Ndn0Y5wWCL9Q^&^k^jM{S3#&(K^Sik>XiufN<>pOzTHr9Cy>UrS-Sd&XCU(S3)HJV$m>3I zVSR_)U&<-9nRV*T-tz?GCR9GFv%tTgH@(pIT-r>_-ZZ*=&Z=>yOR1+~P$$$i~cU_uP|?INxaahG0)xG_cIG=e*Bprz3-E{BL%nuLr&vxi+*L zJ#ljF^_FWq=!PMK6f>+m=vNnBKdeu3r=_VsM4vzW8Ex}2F;DQtXH+z;NuATWu5?Sw zJsoAOx>BuflKb)nC%8%Au;?@X9cdM9=P|r}N4ovV-8FSPSz5lKlFgQ~ zbglH8@2sp?)X5H1`s&@-yy@0pG%C#)VTeco;9 zRM{^-=l{@#t}~2ot!&-~^V*tvbiN$)x>0L7{m3f2l5bnlzNy38jcM76u6FEjvw4Xt zefQ(Fvfun%>9EALTKVrT^v8yysq-Kg>a6p!{`E>rnwjwUuq?ABO>MF>{7I9R^!M#w zHCeu*1^p>gJiD8y1#RilJoZ#|bK1E2YPD3=9PwyQ3nOQYm+Wgs=UEJX{vfRx?YQ5w zYhv4GsO!!&Jm%u`py|%^yvx>Yn|nDUK27QS)yXXuuWm}c#)iGTqiBjaHKh-o-aF`S zI$^y#Azn`O(&TTg54$_jmTO<=18bX5b@Z=AcQTvM>gQ7=3*S1@&sT=H_sMspGr}7h zZnbu#*B>3)Ua-C~4J?{FrFTGM#L?j<(Yla=2__cDnFSVmJpZi5!JY-A9Jhq<{ zHN}?hUs@en-`E!MZAc$S(wgDZv zFMoG#MFXs>2K4Jev{Q2H28g#cRqs68dh$qX#M>J2wxZpf5;q(lZ$-bUbW=~TvqBzP zqK;b9FYGp09P4UHk5^Ci5nZsLDT{_mE9P0y`k_Aa-9jy}?ikwkX05?`BZD|IG+0$2 zzTJYM>&NJsibd~)2ZH7S*h+_4dmATBYAyh7|7!ZMYUF448#@eVPq%R9uBH+2~IZDLxFx5?}k9X*>KyhSAb=h`?0-6H-+3K})} z=_WC)!<%Go1KSB6l{ZMU(*3VO18)%1`nyirKIr>uLezCqJkMr);i+pF?=@n3fAP|W z7S~Ar+`+9n@m=(sK1Ba$;m^`49Im64Uy*NPZi!ExKp%1#$Uj43oOg}+ zx!D;qcDd^CcN0%zeovDtFW=vFSbvHnf46hAyYLhVkDq5h>GVlr%Kww3@gYHESN0@X zkpA#s(7F>OtLtxyB?C^7gxWcF_Z}Z7^xKi^0}GFn$z%E`jt)Ccq+J@@xiQB{=#`75 z*(G}7zO*djaHO7$9-VJH)KO2CMO2E1ryV05?v2%FKddIVM!7h#y{gG*i?j1c)lu^B z`z1|=gdfGcA3;2hkcqQT&s~4%$=hGG!fMN2^5r~c`uNd%iA8kr{PQjLlA8@x zKHk6VA=SsqUDm4ikR2t7p@DP{iBT&HcfH(A{5=DEwp_WJR9(z`8ym2jyjgiUd*_|s z$<4I(BYUj+olNL1&{dB5os71>d~uNb@1()ak=f7dcadck83%K=?IP30Z7hG3zKi6H zJnHsC#4hC7E>iA%bnw&LRpjma{MGx8RgtNy2e@5}uOcZ6E896A*h%(XJ^MWI+HYiP zufp|vY=0wke*6dz?;S+^$F{mP(qEBhzmO9G%t?+8ue|ZJT2zszM zwQ~h2I`2TotkjXLZ_*bp_tlZ7gA%%_HRYsyz00LDOE)1;Hj*wwmjoBzDkBSbFD~`% zQbwZN3?5$C{|EBR^tfFfvo??mmW$5csa;PZG8T#^c34m3cTWtZ*`*{e%~|v2$U0Jb zV^IFDoz{^9()lG5$FC)~FYAA!$4gKLN=W{@V?$ia*N{*Sw?zZUlinMwT!HZTC}b>=v$I~shM`JZYjAD-t0X4sE7>O>@3>i zRfN1PB=gJWtNku4A<499V8_TMsAr4G>K+AOuBu-|vi3FkVqWSZtm}oO)h)-fqeB;x z<$2Zd-yB{*&c_Emdls>PxOCaNP_}PA@@qahxA<0t;ZOncser_cl66_PcOJRZfAVer b+jEgmbBW~jn*s0r=8(W0NrqmD`Q%>!4ooe^ literal 0 HcmV?d00001 diff --git a/data/cross_sections/CLO_JPL06.nc b/data/cross_sections/CLO_JPL06.nc new file mode 100644 index 0000000000000000000000000000000000000000..9b0c4281a9b264c284a894433907027ef90a39e9 GIT binary patch literal 9189 zcmeHNZE#dq89uw4gk)Jh6xx7PIf6w(cW-w$A0(NJY}jl*$S2u+Q9CX-yEoa}+`Tuw z_iku9Fa+yZaL{3x($YU2#A2sb+8C3Fr)2{R>4+$-+RtU zLa16gt>ce9$?miFe7@&-pYxubb9Qe_W77ivXZ+sMQX1t{0 z@dkeAUFGKaig}S}G8g?`;=wAn+(wE3d?jRV4e=bPzpjNRv7oQ)8_K{P5*6;{#FJ@T$?`L=iR`rn*Qi^@Y+Z&QevIMel zDrjqCYe%9p+1k~aO2!-76C`jKPRK`!S%c1oj)Yr3#_RWHHJie}LJif4H#X6XmQxK| zGYvYFcc@`HGezTyr71(YDp8&)Ye)bL@52w^Vx*&KT0$yLP|w$$KwM_U)V3ccb@2 zAQZ^GM>oP!#}p$~^GVe;kbm()T%NvlEFa4rYcBmXqR&GQlV@LYeb7dB0+-&dku$dj=^*iP9P8OeZopJZg z3QtTzplbjZ9+Eiw8O=HW!V|xpEkqHCg{c9Z32_k#qx5?cO^KdFXFQQ|rC$TXpwKdk z^O1!p<}4_?RM-6A2U)bAR|xDtVdG-_xGWM}37oq$goQ#Fm$n?aj24#Qiih|wEH4_; z4EwLb%rIL=^VpZFhP;zCE%)0L(CBKdx38Zz0&#S3*tBx0MI%x;0IhJ{m$FuAONCgO zahu_;iWWx0(W+255{iUrq`E#_U0+=rsAQDzHS1`OuDWfO@iuYP8|8YVGkOcG!>~F1 zrF2;1?a-$<)?hLfZ|_QN?rsm)NX~uEFE)4F3E748|H(jeysO<6Df$nR5UkkKg*JuI zMZWl^DYnA7`9k3$U-ZL9A!Pp|@u`3X1_}%m7$`7MV4%Q2fq?=81qKQX6c{KlP+*|I z!2g7SDI>$NXA$>pB)gc-Kr{9MSU5ZxP_wOaLrnL ztMK>ETrAa>MeAzY`~%aM71UCU5$HNvRTHX=)EK$7U|EHKlQyJl=7^2 zG1H>RJw=5>RN7~&R5NIUlQZqHtO`Fkph+j6fqQI;CdL&#ui$|@?PM;CRKIv2NP%kF zod>egtT0>?ibAgHP`m}t>+ThNyE0_L#hz(Rz_b-x0kU8!=p4x(H*;`P@7TbyaGNU* z%mt;ufKbA{LycU}?+4#m%7|)H6F$&r$jE`!aODdYLf$ZSvMNoRhNBsIGtVYGgs#g_ zdkpRjtbC5`KO-WbF5W9Zs#|6 z$~xquZL2r+)OE|Tl?(2kd~1uW9~o^k4kYDg79Ib7$+A9q@JeLfLuSA1y)ysAnf(Lu zq59)LIHe5Af4Mf5`OO=H@|8CqeCxHh2jzF7`@i+@#8&xKeY$=3i2?bQ9s9m~aQT3I zBK7XkbKCpn=dSFo_1)elA2Ze*-J$i#zux&z`FGDIWAcU4yIxJ*8h8DQ(1rWar^Z@+#T>t<8 literal 0 HcmV?d00001 diff --git a/data/cross_sections/HNO3_JPL06.nc b/data/cross_sections/HNO3_JPL06.nc new file mode 100644 index 0000000000000000000000000000000000000000..f56b8dd432d353fcbcb9d2219326f93da1c8f796 GIT binary patch literal 3436 zcmb7_dsGuw9>=5NgD6@RySUZ$;;TWEGzmo>H};O8f*`N(7J&@OKt_^Gm`o6~)+iRO zcD1W})LI|7QdhT%XnpPx4TpoaKCnIjAJs*xwN~qPeNkQ1y)(J7cy{-X=A6lAe!qL? zcYpW$`^}x0$y2+@ygxZc4}AF*vUVqC3clnHnx}1y!0@6+ChLx`egb22Fgz`|c*g6+ zKJ53ySFcGf!NT$Q&O0|SJtGrLrUe>g7je9e;el4I!DGEmJVOhN**9;nPNREQrO~Ri z8lVk}(u75YMTlcP1y-;!;+;2MPz=vGoxsVM1eUXd*Rt~~Vx_5KjsvsP3gVKILkr?l zrfVjsg>pd(A?r)x$Y;}KjFqt$3l?d>r-8T+@xG_a&I)96V!ZBln>4Q%zIEf_d)Nvx ziMbSzTnb*xMSP}fOwsOPF-yYWd%q-IIEc^y-)|HsVtJmnl(aZZPY5yedtTxOPWs=<^p+UlaER!s6%OB>d>QVFos zAjO5nGI<9z6I43vWdxAIa}Jl)DK_tPx{$+f5yw~HZPQK~i-or$7Dp|Q;cPa>E;zBu zzOqSuA$oC>;K%U@sSdE)lz=XDa=b&zL~2XlE8+x5+FX}+7MM7usE9SO;sPKVHyB<- z9jyk@F22Mf7ELNANXL6CG1J)iAC(XVOT*RYtelO*qK@p+PQc7KaGH$W#E9K1SQsF_ zivS-5|2IwLNPtTi%W$~0u&@SE-N*)^{TcMej-dHDFY@>m~WXsou$blHw-~YD<#HY)DpILrlQ0&KN zR!@vRR?2(U$tVkjZ5gJU(Bma5K}N9@vE76oFK0=RQL`vwy9qsBmPn9MB^0sUgdQ)g z5@eK>BDR~*R7nOC`vtQb25q&`9VetS0mjwtMOQF3Kn#5St=261oYi2|a}EUXl+I)BS;B=*c73W&`nrP=pk(P(j`Ggxd5>#LL;G@u$s_A*e+3TMI*LP z>)51cQGK|!C{oYtINrK6dcU4-h#mUy*e?B)Mxt8ccL#(XsOpLZfmbZPioKuo&MUUs`2_E%MRt1Z+}}KQ2BN8 z{&PF@J)aL)wNG|Z-;KTVbjoZM?93aMHT3=&*zraE`gMut;iKTbmk)fh7ykUqss{fN zAHjP^GhQY|HNqbpnxDGgIt;&SnI~V784Ir`)LjZIXolBZD^tS%xCLIR4o})Q^9a1; z4h>j!a1gwh+HvYZa~^D4zq#Ma=uxn3(&HU;XK&ax>hQ;@H@CsIo|V=~r>DJQ>(>L7 z!5f#s){i&d>&Bgit%Ltw(lE#XTh_JQb?L zb8n0e+0bzSo|TRC4;VwkGsU?n2jf12r{y=#l)1Fq+s^fXwx)sZc_qtX@s!oz+R4q(Xy~Wi8j%OHzbZ}t$CKmG zP<3UEs&6Zd4=DQnvGr#d6E{t(+H8Rnjqm5=GcFju`*tUcn*p`scKuMGZH9{5<}od2 zC&SSvM>FbE%OR+}S+n?f8w^r?Ne$YQ3ImgKg0l)&`zSSGmp_zyIdy1r?*f$DXI*#K zn5igNx6D>`F9_w<71SQ4mZO;q-#j*C#?p zTIMNTTIW?MrBas`e1RCA<1axLgJ`szDxq1g-1FDiSm70oG0 zse3T@OEkCCwLExr8k(1O`l{W}gFd)-ba!qRLT-|0Wk>U+Hy4jbi?FIMswz?C_HmjF u9gRLD@v2CkOT=XYenLK0_HDm*Zb>s*a{R_Z>VO+9Dd_d&&r=0N{Qm)L&@qqz literal 0 HcmV?d00001 diff --git a/examples/ts1_tsmlt.json b/examples/ts1_tsmlt.json index b00663f3..fb668cc0 100644 --- a/examples/ts1_tsmlt.json +++ b/examples/ts1_tsmlt.json @@ -321,7 +321,7 @@ "__reaction": "HNO3 + hv -> OH + NO2", "cross section": { "netcdf files": [ - { "file path": "data/cross_sections/HNO3_1.nc" } + { "file path": "data/cross_sections/HNO3_JPL06.nc" } ], "type": "HNO3+hv->OH+NO2" }, @@ -555,11 +555,7 @@ "cross section": { "netcdf files": [ { - "file path": "data/cross_sections/BrO_1.nc", - "interpolator": { - "type": "fractional target", - "fold in": true - } + "file path": "data/cross_sections/BRO_JPL06.nc" } ], "type": "base" @@ -818,7 +814,7 @@ "__comments": "TODO - this doesn't exactly match the products in TS1", "cross section": { "netcdf files": [ - { "file path": "data/cross_sections/ClOOCl_1.nc" } + { "file path": "data/cross_sections/CL2O2_JPL10.nc" } ], "type": "base" }, @@ -840,6 +836,20 @@ "type": "ClO+hv->Cl+O(1D)" } }, + { + "name": "jclo_o3p", + "__reaction": "ClO + hv -> Cl + O", + "cross section": { + "netcdf files": [ + { "file path": "data/cross_sections/CLO_JPL06.nc" } + ], + "type": "base" + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + } + }, { "name": "jclono2_a", "__reaction": "ClONO2 + hv -> Cl + NO3", diff --git a/test/data/xsqy.doug.config.json b/test/data/xsqy.doug.config.json index 901c9c2d..971775d8 100644 --- a/test/data/xsqy.doug.config.json +++ b/test/data/xsqy.doug.config.json @@ -109,5 +109,126 @@ }, "label": "CH2Br2 + hv -> 2Br", "tolerance": 1.0e-3 + }, + { + "cross section": { + "type": "base", + "netcdf files": [ + { "file path": "data/cross_sections/BRO_JPL06.nc" } + ] + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + }, + "label": "BRO + hv -> Br + O", + "__note": "first test: excluding edges of interpolation because of double vs float algoritms", + "mask" : [ { "index": 62 }, { "index": 86 }] + }, + { + "cross section": { + "type": "base", + "netcdf files": [ + { "file path": "data/cross_sections/BRO_JPL06.nc" } + ] + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + }, + "label": "BRO + hv -> Br + O", + "__note": "second test: including edges of interpolation with relaxed tolerance", + "tolerance": 5.0e-3 + }, + { + "cross section": { + "type": "base", + "netcdf files": [ + { "file path": "data/cross_sections/CL2O2_JPL10.nc" } + ] + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + }, + "label": "Cl2O2 + hv -> Cl + ClOO", + "__note": "first test: excluding edges of interpolation because of double vs float algoritms", + "mask" : [ { "index": 34 }, { "index": 97 } ] + }, + { + "cross section": { + "type": "base", + "netcdf files": [ + { "file path": "data/cross_sections/CL2O2_JPL10.nc" } + ] + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + }, + "label": "Cl2O2 + hv -> Cl + ClOO", + "__note": "second test: including edges of interpolation with relaxed tolerance", + "tolerance": 1.0e-3 + }, + { + "cross section": { + "type": "base", + "netcdf files": [ + { "file path": "data/cross_sections/CLO_JPL06.nc" } + ] + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + }, + "label": "ClO + hv -> Cl + O", + "__note": "first test: excluding edges of interpolation because of double vs float algoritms", + "mask": [ { "index": 51 }, { "index": 71 }] + }, + { + "cross section": { + "type": "base", + "netcdf files": [ + { "file path": "data/cross_sections/CLO_JPL06.nc" } + ] + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + }, + "label": "ClO + hv -> Cl + O", + "__note": "second test: including edges of interpolation with relaxed tolerance", + "tolerance": 1.0e-3 + }, + { + "cross section": { + "type": "HNO3+hv->OH+NO2", + "netcdf files": [ + { "file path": "data/cross_sections/HNO3_JPL06.nc" } + ] + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + }, + "label": "HNO3 + hv -> OH + NO2", + "__note": "first test: excluding edges of interpolation because of double vs float algoritms", + "mask": [ { "index": 30 }, { "index": 79 } ] + }, + { + "cross section": { + "type": "HNO3+hv->OH+NO2", + "netcdf files": [ + { "file path": "data/cross_sections/HNO3_JPL06.nc" } + ] + }, + "quantum yield": { + "type": "base", + "constant value": 1.0 + }, + "label": "HNO3 + hv -> OH + NO2", + "__note": "second test: including lower edge of interpolation with relaxed tolerance (upper edge is a very small value with large relative difference)", + "tolerance": 1.0e-3, + "mask": [ { "index": 79 } ] } ] diff --git a/test/unit/tuv_doug/JCALC/CMakeLists.txt b/test/unit/tuv_doug/JCALC/CMakeLists.txt index 93fd2604..6e70ac66 100644 --- a/test/unit/tuv_doug/JCALC/CMakeLists.txt +++ b/test/unit/tuv_doug/JCALC/CMakeLists.txt @@ -3,8 +3,12 @@ target_sources(tuv_doug PRIVATE + XSQY_BRO.f XSQY_CH2BR2.f + XSQY_CL2O2.f + XSQY_CLO.f XSQY_H2O.f + XSQY_HNO3.f XSQY_N2O5.f ) diff --git a/test/unit/tuv_doug/JCALC/XSQY_BRO.f b/test/unit/tuv_doug/JCALC/XSQY_BRO.f new file mode 100644 index 00000000..df054c6e --- /dev/null +++ b/test/unit/tuv_doug/JCALC/XSQY_BRO.f @@ -0,0 +1,116 @@ + subroutine XSQY_BRO(nw,wl,wc,nz,tlev,airlev,j,sq,jlabel,pn) +!-----------------------------------------------------------------------------! +! purpose: ! +! provide product (cross section) x (quantum yield): ! +! BrO + hv -> Br + O ! +! cross section: JPL06 ! +! quantum yield: is unity. ! +!-----------------------------------------------------------------------------! +! parameters: ! +! nw - integer, number of specified intervals + 1 in working (i) ! +! wavelength grid ! +! wl - real, vector of lower limits of wavelength intervals in (i) ! +! working wavelength grid ! +! wc - real, vector of center points of wavelength intervals in (i) ! +! working wavelength grid ! +! nz - integer, number of altitude levels in working altitude grid (i) ! +! tlev - real, temperature (k) at each specified altitude level (i) ! +! airlev - real, air density (molec/cc) at each altitude level (i) ! +! j - integer, counter for number of weighting functions defined (io) ! +! sq - real, cross section x quantum yield (cm^2) for each (o) ! +! photolysis reaction defined, at each defined wavelength and ! +! at each defined altitude level ! +! jlabel - character*60, string identifier for each photolysis reaction (o) ! +! defined ! +!-----------------------------------------------------------------------------! +! edit history: ! +! 07/27/07 Doug Kinnison ! +!-----------------------------------------------------------------------------! + implicit none + include 'params' + +!-----------------------------------------------------------------------------! +! ... input ! +!-----------------------------------------------------------------------------! + real, intent(in) :: wl(kw) + real, intent(in) :: wc(kw) + real, intent(in) :: tlev(kz) + real, intent(in) :: airlev(kz) + + integer, intent(in) :: nz + integer, intent(in) :: nw + + character*80, intent(in) :: pn + character*60, intent(out) :: jlabel(kj) + real, intent(out) :: sq(kj,kz,kw) + +!-----------------------------------------------------------------------------! +! ... input/output ! +!-----------------------------------------------------------------------------! + integer, intent(inout) :: j + +!-----------------------------------------------------------------------------! +! ... local ! +!-----------------------------------------------------------------------------! + integer kdata + parameter(kdata=300) + integer i, iw, n, idum, ierr, iz + real x1(kdata) + real y1(kdata) + real yg(kw) + real qy + +!---------------------------------------------- +! ... jlabel(j) = 'BRO + hv -> Br + O' +!---------------------------------------------- + j = j+1 + jlabel(j) = 'BRO + hv -> Br + O' + +!---------------------------------------------------- +! ... cross sections from JPL06 recommendation +!---------------------------------------------------- +! ... 0.5nm resolution JPL06. + open(kin,file=TRIM(pn)//'XS_BRO_JPL06.txt',status='old') + + read(kin,*) idum, n + do i = 1, idum-2 + read(kin,*) + enddo + + do i = 1, n + read(kin,*) x1(i), y1(i) + enddo + close(kin) + + call addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.) + call addpnt(x1,y1,kdata,n, 0.,0.) + call addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.) + call addpnt(x1,y1,kdata,n, 1e38,0.) + call inter2(nw,wl,yg, n,x1, y1,ierr) + + if (ierr .ne. 0) then + write(*,*) ierr, jlabel(j) + stop + endif + +!------------------------------------------------------- +! ... quantum yield (assumed) to be unity (JPL06) +!------------------------------------------------------- + qy = 1.0 + + do iw = 1, nw-1 + do iz = 1, nz + sq(j,iz,iw) = qy * yg(iw) + + enddo + enddo + +!------------------------------------------------------- +! ... Check Routine +! do iw = 61, 87 +! print*, iw, wc(iw), (qy * yg(iw)) +! enddo +! stop +!------------------------------------------------------- + + end subroutine XSQY_BRO diff --git a/test/unit/tuv_doug/JCALC/XSQY_CL2O2.f b/test/unit/tuv_doug/JCALC/XSQY_CL2O2.f new file mode 100644 index 00000000..87abe738 --- /dev/null +++ b/test/unit/tuv_doug/JCALC/XSQY_CL2O2.f @@ -0,0 +1,117 @@ + subroutine XSQY_CL2O2(nw,wl,wc,nz,tlev,airlev,j,sq,jlabel,pn) +!-----------------------------------------------------------------------------! +! purpose: ! +! provide product (cross section) x (quantum yield) for hcl photolysis: ! +! cl2o2 + hv -> cl + cloo ! +! cross section: from JPL10 ! +! ! +!-----------------------------------------------------------------------------! +! parameters: ! +! nw - integer, number of specified intervals + 1 in working (i) ! +! wavelength grid ! +! wl - real, vector of lower limits of wavelength intervals in (i) ! +! working wavelength grid ! +! wc - real, vector of center points of wavelength intervals in (i) ! +! working wavelength grid ! +! nz - integer, number of altitude levels in working altitude grid (i) ! +! tlev - real, temperature (k) at each specified altitude level (i) ! +! airlev - real, air density (molec/cc) at each altitude level (i) ! +! j - integer, counter for number of weighting functions defined (io) ! +! sq - real, cross section x quantum yield (cm^2) for each (o) ! +! photolysis reaction defined, at each defined wavelength and ! +! at each defined altitude level ! +! jlabel - character*60, string identifier for each photolysis reaction (o) ! +! defined ! +!-----------------------------------------------------------------------------! +! edit history: ! +! 01/13/2012 Doug Kinnison ! +!-----------------------------------------------------------------------------! + implicit none + include 'params' + +!-----------------------------------------------------------------------------! +! ... input ! +!-----------------------------------------------------------------------------! + real, intent(in) :: wl(kw) + real, intent(in) :: wc(kw) + real, intent(in) :: tlev(kz) + real, intent(in) :: airlev(kz) + + integer, intent(in) :: nz + integer, intent(in) :: nw + + character*80, intent(in) :: pn + character*60, intent(out) :: jlabel(kj) + real, intent(out) :: sq(kj,kz,kw) + +!-----------------------------------------------------------------------------! +! ... input/output ! +!-----------------------------------------------------------------------------! + integer, intent(inout) :: j + +!-----------------------------------------------------------------------------! +! ... local ! +!-----------------------------------------------------------------------------! + integer kdata + parameter(kdata=600) + integer i, iw, n, idum, ierr, iz + real x1(kdata) + real y1(kdata) + real yg(kw) + real qy + +!---------------------------------------------- +! ... jlabel(j) = 'cl2o2 -> cl + cloo' +!---------------------------------------------- + j = j+1 + jlabel(j) = 'Cl2O2 + hv -> Cl + ClOO' +! print*,jlabel(j) +!---------------------------------------------------- +! ... cross sections +!---------------------------------------------------- + open(kin,file= + $ TRIM(pn)//'XS_CL2O2_JPL10_500nm.txt',status='old') + + read(kin,*) idum, n + do i = 1, idum-2 + read(kin,*) + enddo + do i = 1, n + read(kin,*) x1(i), y1(i) + enddo + close(kin) + +!---------------------------------------------------------- +! do i = 1, n +! print*, i, x1(i), y1(i) +! enddo +! stop +!---------------------------------------------------------- + call addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.) + call addpnt(x1,y1,kdata,n, 0.,0.) + call addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.) + call addpnt(x1,y1,kdata,n, 1e38,0.) + call inter2(nw,wl,yg,n,x1,y1,ierr) + + + if (ierr .ne. 0) then + write(*,*) ierr, jlabel(j) + stop + endif + +!-------------------------------------------------------------- +! ... quantum yield assumed to be 1.0 +!-------------------------------------------------------------- + qy = 1.0 + do iw = 1, nw-1 + do iz = 1, nz + sq(j,iz,iw) = qy * yg(iw) + enddo + enddo +!---------------------------------------------------------- +! do iw = 28, 99 +! print*, iw, wc(iw), yg(iw) +! enddo +! stop +!---------------------------------------------------------- + end subroutine XSQY_CL2O2 diff --git a/test/unit/tuv_doug/JCALC/XSQY_CLO.f b/test/unit/tuv_doug/JCALC/XSQY_CLO.f new file mode 100644 index 00000000..2dab1cae --- /dev/null +++ b/test/unit/tuv_doug/JCALC/XSQY_CLO.f @@ -0,0 +1,115 @@ + subroutine XSQY_CLO(nw,wl,wc,nz,tlev,airlev,j,sq,jlabel,pn) +!-----------------------------------------------------------------------------! +! purpose: ! +! provide product (cross section) x (quantum yield): ! +! ClO + hv -> Cl + O ! +! cross section: JPL06 ! +! quantum yield: is unity. ! +!-----------------------------------------------------------------------------! +! parameters: ! +! nw - integer, number of specified intervals + 1 in working (i) ! +! wavelength grid ! +! wl - real, vector of lower limits of wavelength intervals in (i) ! +! working wavelength grid ! +! wc - real, vector of center points of wavelength intervals in (i) ! +! working wavelength grid ! +! nz - integer, number of altitude levels in working altitude grid (i) ! +! tlev - real, temperature (k) at each specified altitude level (i) ! +! airlev - real, air density (molec/cc) at each altitude level (i) ! +! j - integer, counter for number of weighting functions defined (io) ! +! sq - real, cross section x quantum yield (cm^2) for each (o) ! +! photolysis reaction defined, at each defined wavelength and ! +! at each defined altitude level ! +! jlabel - character*60, string identifier for each photolysis reaction (o) ! +! defined ! +!-----------------------------------------------------------------------------! +! edit history: ! +! 07/27/07 Doug Kinnison ! +!-----------------------------------------------------------------------------! + implicit none + include 'params' + +!-----------------------------------------------------------------------------! +! ... input ! +!-----------------------------------------------------------------------------! + real, intent(in) :: wl(kw) + real, intent(in) :: wc(kw) + real, intent(in) :: tlev(kz) + real, intent(in) :: airlev(kz) + + integer, intent(in) :: nz + integer, intent(in) :: nw + + character*80, intent(in) :: pn + character*60, intent(out) :: jlabel(kj) + real, intent(out) :: sq(kj,kz,kw) + +!-----------------------------------------------------------------------------! +! ... input/output ! +!-----------------------------------------------------------------------------! + integer, intent(inout) :: j + +!-----------------------------------------------------------------------------! +! ... local ! +!-----------------------------------------------------------------------------! + integer kdata + parameter(kdata=300) + integer i, iw, n, idum, ierr, iz + real x1(kdata) + real y1(kdata) + real yg(kw) + real qy + +!---------------------------------------------- +! ... jlabel(j) = 'ClO + hv -> Cl + O' +!---------------------------------------------- + j = j+1 + jlabel(j) = 'ClO + hv -> Cl + O' + +!---------------------------------------------------- +! ... cross sections from JPL06 recommendation +!---------------------------------------------------- + open(kin,file=TRIM(pn)//'XS_CLO_JPL06.txt',status='old') + + read(kin,*) idum, n + do i = 1, idum-2 + read(kin,*) + enddo + + do i = 1, n + read(kin,*) x1(i), y1(i) + enddo + close(kin) + + call addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.) + call addpnt(x1,y1,kdata,n, 0.,0.) + call addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.) + call addpnt(x1,y1,kdata,n, 1e38,0.) + + call inter2(nw,wl,yg,n,x1,y1,ierr) + + if (ierr .ne. 0) then + write(*,*) ierr, jlabel(j) + stop + endif +!------------------------------------------------------- +! ... quantum yield (assumed) to be unity (JPL06) +!------------------------------------------------------- + qy = 1.0 + + do iw = 1, nw-1 + do iz = 1, nz + sq(j,iz,iw) = qy * yg(iw) + + enddo + enddo +!------------------------------------------------------- +! ... Check routine (no temperature dependence +! print*,'jclo' +! do iw = 30, 72 +! print*, iw, wc(iw), (qy * yg(iw)) +! enddo +! stop +!------------------------------------------------------- + + end subroutine XSQY_CLO diff --git a/test/unit/tuv_doug/JCALC/XSQY_HNO3.f b/test/unit/tuv_doug/JCALC/XSQY_HNO3.f new file mode 100644 index 00000000..67792b40 --- /dev/null +++ b/test/unit/tuv_doug/JCALC/XSQY_HNO3.f @@ -0,0 +1,140 @@ + subroutine XSQY_HNO3(nw,wl,wc,nz,tlev,airlev,j,sq,jlabel,pn) +!-----------------------------------------------------------------------------! +! purpose: ! +! provide product of (cross section) x (quantum yield) for photolysis ! +! hno3 + hv -> oh + no2 ! +! cross section: burkholder et al., 1993 (and JPL06) ! +! quantum yield: assumed to be unity ! +!-----------------------------------------------------------------------------! +! parameters: ! +! nw - integer, number of specified intervals + 1 in working (i) ! +! wavelength grid ! +! wl - real, vector of lower limits of wavelength intervals in (i) ! +! working wavelength grid ! +! wc - real, vector of center points of wavelength intervals in (i) ! +! working wavelength grid ! +! nz - integer, number of altitude levels in working altitude grid (i) ! +! tlev - real, temperature (k) at each specified altitude level (i) ! +! airlev - real, air density (molec/cc) at each altitude level (i) ! +! j - integer, counter for number of weighting functions defined (io) ! +! sq - real, cross section x quantum yield (cm^2) for each (o) ! +! photolysis reaction defined, at each defined wavelength and ! +! at each defined altitude level ! +! jlabel - character*60, string identifier for each photolysis reaction (o) ! +! defined ! +!-----------------------------------------------------------------------------! +! edit history: ! +! 05/98 original, adapted from former jspec1 subroutine ! +! 01/15/08 minor update,dek ! +!-----------------------------------------------------------------------------! + implicit none + include 'params' + +!-----------------------------------------------------------------------------! +! ... input ! +!-----------------------------------------------------------------------------! + real, intent(in) :: wl(kw) + real, intent(in) :: wc(kw) + real, intent(in) :: tlev(kz) + real, intent(in) :: airlev(kz) + + integer, intent(in) :: nz + integer, intent(in) :: nw + + character*80, intent(in) :: pn + character*60, intent(out) :: jlabel(kj) + real, intent(out) :: sq(kj,kz,kw) + +!-----------------------------------------------------------------------------! +! ... input/output ! +!-----------------------------------------------------------------------------! + integer, intent(inout) :: j + +!-----------------------------------------------------------------------------! +! ... local ! +!-----------------------------------------------------------------------------! + integer kdata + parameter(kdata=100) + integer n1, n2 + integer i, iw, n, idum, iz + integer ierr + real x1 (kdata), x2 (kdata) + real y1 (kdata), y2 (kdata) + real yg1(kw), yg2(kw) + real yg( kw) + real tin(nz) + +!---------------------------------------------- +! ... tin set to tlev +!---------------------------------------------- + tin(:) = tlev(:) + +!---------------------------------------------- +! ... jlabel(j) = 'HNO3 -> OH + NO2 +!---------------------------------------------- + j = j + 1 + jlabel(j) = 'HNO3 + hv -> OH + NO2' + +!----------------------------------------------------------------------- +! ... hno3 cross section parameters from burkholder et al. 1993 +!----------------------------------------------------------------------- + open(kin,file=TRIM(pn)//'XS_HNO3_JPL06.txt',status='old') + +!... read lambda and cross sections + read(kin,*) idum, n + do i = 1, idum-2 + read(kin,*) + enddo + do i = 1, n + read(kin,*) x1(i), y1(i) + enddo + +!... read lambda and T-dep coeff. + read(kin,*) + do i = 1, n + read(kin,*) x2(i), y2(i) + enddo + close(kin) + + call addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.) + call addpnt(x1,y1,kdata,n, 0.,0.) + call addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.) + call addpnt(x1,y1,kdata,n, 1e38,0.) + call inter2(nw,wl,yg1,n,x1,y1,ierr) + if (ierr .ne. 0) then + write(*,*) ierr, jlabel(j) + stop + endif + + n= 80 + call addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.) + call addpnt(x2,y2,kdata,n, 0.,0.) + call addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),0.) + call addpnt(x2,y2,kdata,n, 1.e+38,0.) + call inter2(nw,wl,yg2,n,x2,y2,ierr) + if (ierr .ne. 0) then + write(*,*) ierr, jlabel(j) + stop + endif + +!-------------------------------------------------- +! ... quantum yield = 1 +! correct for temperature dependence +!-------------------------------------------------- + do iw = 1, nw - 1 + do iz = 1, nz + sq(j,iz,iw) = yg1(iw) + $ * exp( yg2(iw)/1.e3*(tin(iz)-298.) ) + enddo + enddo + +!------------------------------------------------------- +! ... Check routine (no temperature dependence +! iz = 1 +! do iw = 29, 79 +! print*, iw, wc(iw), sq(j,iz,iw) +! enddo +! stop +!------------------------------------------------------- + + end subroutine XSQY_HNO3 diff --git a/test/unit/tuv_doug/data_sets.F90 b/test/unit/tuv_doug/data_sets.F90 index 434014ee..50874616 100644 --- a/test/unit/tuv_doug/data_sets.F90 +++ b/test/unit/tuv_doug/data_sets.F90 @@ -45,9 +45,10 @@ subroutine test_data_set( ) class(cross_section_t), pointer :: cross_section class(quantum_yield_t), pointer :: quantum_yield - character(len=*), parameter :: Iam = "H2O cross section test" + character(len=*), parameter :: Iam = "Doug's cross section tests" type(config_t) :: config, config_pair, cs_config, qy_config - class(iterator_t), pointer :: iter + type(config_t) :: mask_points_config, mask_point_config + class(iterator_t), pointer :: iter, mask_points_iter type(string_t) :: cs_type_name, qy_type_name, label character, allocatable :: buffer(:) integer :: pos, pack_size @@ -59,7 +60,9 @@ subroutine test_data_set( ) class(profile_t), pointer :: air, temperature class(grid_t), pointer :: wavelength real(kind=dk) :: tolerance + integer, allocatable :: mask_points(:) integer :: i + logical :: found ! Load grids based on Doug's TUV grids => get_grids( ) @@ -84,6 +87,21 @@ subroutine test_data_set( ) call config_pair%get( "label", label, Iam ) call config_pair%get( "tolerance", tolerance, Iam, & default = 1.0e-6_dk ) + call config_pair%get( "mask", mask_points_config, Iam, & + found = found ) + if (found) then + mask_points_iter => mask_points_config%get_iterator( ) + allocate( mask_points( mask_points_config%number_of_children( ) ) ) + do i = 1, size( mask_points ) + call assert( 564855121, mask_points_iter%next( ) ) + call mask_points_config%get( mask_points_iter, mask_point_config, & + Iam ) + call mask_point_config%get( "index", mask_points( i ), Iam ) + end do + call assert( 888375064, .not. mask_points_iter%next( ) ) + else + allocate(mask_points(0)) + end if ! Load and test cross section if( musica_mpi_rank( comm ) == 0 ) then @@ -144,6 +162,8 @@ subroutine test_data_set( ) ! Skip first two bins because Lyman-Alpha bins are different in ! Doug's version of TUV-x. Data sets were adapted to have Lyman-Alpha ! specific data go into the TUV-x Lyman-Alpha bin 121.4-121.9 nm + ! Also skip any points explicitly masked in the configuration + tuvx_xsqy(:,mask_points(:)) = doug_xsqy(:,mask_points(:)) call check_values( 377150482, tuvx_xsqy(:,3:), & real( doug_xsqy(:,3:), kind=dk ), tolerance ) @@ -152,6 +172,7 @@ subroutine test_data_set( ) deallocate( cross_section_data ) deallocate( quantum_yield_data ) deallocate( tuvx_xsqy ) + deallocate( mask_points ) end do diff --git a/test/unit/tuv_doug/driver.F90 b/test/unit/tuv_doug/driver.F90 index cbdd8a7c..1dea6336 100644 --- a/test/unit/tuv_doug/driver.F90 +++ b/test/unit/tuv_doug/driver.F90 @@ -141,6 +141,18 @@ subroutine calculate( label, temperature, air_density, xsqy ) case( "CH2Br2 + hv -> 2Br" ) call XSQY_CH2BR2(nw,wl,wc,nz,temperature,air_density,j,l_xsqy,all_labels,pn) xsqy(:,:) = l_xsqy(1,:nz,:nw) + case( "BRO + hv -> Br + O" ) + call XSQY_BRO(nw,wl,wc,nz,temperature,air_density,j,l_xsqy,all_labels,pn) + xsqy(:,:) = l_xsqy(1,:nz,:nw) + case( "Cl2O2 + hv -> Cl + ClOO" ) + call XSQY_CL2O2(nw,wl,wc,nz,temperature,air_density,j,l_xsqy,all_labels,pn) + xsqy(:,:) = l_xsqy(1,:nz,:nw) + case( "ClO + hv -> Cl + O" ) + call XSQY_CLO(nw,wl,wc,nz,temperature,air_density,j,l_xsqy,all_labels,pn) + xsqy(:,:) = l_xsqy(1,:nz,:nw) + case ("HNO3 + hv -> OH + NO2" ) + call XSQY_HNO3(nw,wl,wc,nz,temperature,air_density,j,l_xsqy,all_labels,pn) + xsqy(:,:) = l_xsqy(1,:nz,:nw) case default call die( 946669022 ) end select diff --git a/tool/data_conversion/photo.config.json b/tool/data_conversion/photo.config.json new file mode 100644 index 00000000..5ee4685d --- /dev/null +++ b/tool/data_conversion/photo.config.json @@ -0,0 +1,44 @@ +{ + "photoreactions": [ + { + "molecule": "BRO", + "cross-sections": [ + { + "filespec": "XS_BRO_JPL06.txt", + "nPreSkip": 17, + "nRead": 198 + } + ] + }, + { + "molecule": "CL2O2", + "cross-sections": [ + { + "filespec": "XS_CL2O2_JPL10_500nm.txt", + "nPreSkip": 32, + "nRead": 521 + } + ] + }, + { + "molecule": "CLO", + "cross-sections": [ + { + "filespec": "XS_CLO_JPL06.txt", + "nPreSkip": 19, + "nRead": 36 + } + ] + }, + { + "molecule": "HNO3", + "cross-sections": [ + { + "filespec": "XS_HNO3_JPL06.txt", + "nPreSkip": 26, + "nRead": 80 + } + ] + } + ] +} diff --git a/tool/data_conversion/text_to_netcdf.py b/tool/data_conversion/text_to_netcdf.py new file mode 100644 index 00000000..bdf99226 --- /dev/null +++ b/tool/data_conversion/text_to_netcdf.py @@ -0,0 +1,61 @@ +#!/Library/Frameworks/Python.framework/Versions/3.8/bin/python3 + +import numpy as np +import sys +import json +from xsqy_subs import xform_to_netCDF + +#----------------------------------------------------- +# json config file is only possible argument +#----------------------------------------------------- +if( len(sys.argv) > 2 ): + print(f'\n{sys.argv[0]}: requires one or no arguments') + sys.exit( -1) +elif( len(sys.argv) == 2 ): + filespec = sys.argv[1] +else: + filespec = 'photo.config.json' + + +#----------------------------------------------------- +# open json photo config file +#----------------------------------------------------- +#ilespec = 'photo.config.tst.json' +try: + fp = open(filespec,'r') +except: + print(f"Failed to open {filespec}") + sys.exit(-1) + +#----------------------------------------------------- +# transfer config file into dictionary +#----------------------------------------------------- +try: + photDict = json.load(fp) +except: + print(f"Failed to load json file {filespec}") + sys.exit(-1) + +#----------------------------------------------------- +# done with json input file +#----------------------------------------------------- +fp.close() + +#----------------------------------------------------- +# loop through photo reactions in dictionary +#----------------------------------------------------- +list = photDict['photoreactions'] +molecule = "" +for rxt in list: +#----------------------------------------------------- +# call xform_to_netCDF +#----------------------------------------------------- + if( molecule != rxt['molecule']): + nFile = 1 + molecule = rxt['molecule'] + else: + nFile += 1 + xform_to_netCDF(nFile,rxt,'./') + +print('\n') +print(f'\nThere are {len(list)} photoreactions in {filespec}') diff --git a/tool/data_conversion/xsqy_subs.py b/tool/data_conversion/xsqy_subs.py new file mode 100644 index 00000000..d7f53404 --- /dev/null +++ b/tool/data_conversion/xsqy_subs.py @@ -0,0 +1,218 @@ +import numpy as np +import sys +from netCDF4 import Dataset +import netCDF4 as ncd +from datetime import datetime as dt + +""" +Function to read the data file(s) +""" +def read_data_file(data_dictionary,dataTray): + + InpFileSpec = data_dictionary['filespec'] + try: + InpFile = open(InpFileSpec,'r') + except: + print(f'Failed to open data file {InpFileSpec}') + sys.exit(-3) + + print(f'Opened data file {InpFileSpec}') + nLines = len(InpFile.readlines()) + InpFile.seek(0) + nskipHdr = data_dictionary['nPreSkip'] if 'nPreSkip' in data_dictionary else 0 + nRead = data_dictionary['nRead'] if 'nRead' in data_dictionary else 0 + + header = '' +# if header lines exist then read them + if( nskipHdr > 0 ): + for ndx in range(nskipHdr): + header += InpFile.readline() + InpFile.seek(0) + + nskipHdr = abs(nskipHdr) + nskipEnd = nLines - (nskipHdr + nRead) + print(f'nLines,nskipHdr,nRead,nskipEnd = {nLines},{nskipHdr},{nRead},{nskipEnd}') + try: + data = np.genfromtxt(InpFile,dtype='float64',skip_header=nskipHdr,skip_footer=nskipEnd,comments=None) + print(f'Read cross section file {InpFileSpec}') + except: + print(f'Failed to read data file {InpFileSpec}') + sys.exit(-2) + + try: + dataTray.append(data) + except: + print('Failed to append data array to dataTray') + sys.exit(-2) + + InpFile.close() + print(f'Closed data file {InpFileSpec}') + return(header) + +""" +Function to write the netCDF file +""" +def stuff_netCDF_file(ncFile,interpolationTemps,dataTray,hasLambdaGrid,InpFileSpecs,var_typ,Headers): + + ndataVars = len(dataTray) + print(f'ndataVars = {ndataVars}') + + for dataVarNdx in range(ndataVars): + nparameterRow,nparameterCol = np.shape(dataTray[dataVarNdx]) + if( hasLambdaGrid ): + nparameterCol -= 1 + ntemperature = min( len(interpolationTemps),nparameterCol ) + print(f'data array is ({nparameterRow},{nparameterCol})') + DataTag = var_typ + "_parameters" +# define dimensions + RowDimName = 'bins' + ColDimName = 'parameters' + TempDimName = 'temperatures' + print(f'Variable type = {var_typ}') + print(f'RowDimName,ColDimName = {RowDimName},{ColDimName}') + ncFile.createDimension(RowDimName,nparameterRow) + ncFile.createDimension(ColDimName,nparameterCol) + ncFile.createDimension(TempDimName,ntemperature) +# create wavelength grid + if( hasLambdaGrid ): + Var = ncFile.createVariable('wavelength',np.float64,(RowDimName)) + Var.units = 'nm' +# write wavelength grid + Var[:] = dataTray[dataVarNdx][:,0] +# create interpolation temperature array + if( len(interpolationTemps) > 0 ): + Var = ncFile.createVariable('temperature',np.float64,(TempDimName)) + Var.units = 'K' +# write interpolation temperature array + Var[:] = interpolationTemps[:ntemperature] +# create cross section or quantum yield data array + Var = ncFile.createVariable(DataTag,np.float64,(ColDimName,RowDimName)) + Var.hdr = Headers[dataVarNdx] +# write data array + if( hasLambdaGrid ): + Var[:,:] = np.transpose(dataTray[dataVarNdx][:,1:]) + else: + Var[:,:] = np.transpose(dataTray[dataVarNdx]) + if( var_typ == 'cross_section' ): + if( hasLambdaGrid ): + Var.units = 'cm^2' + else: + Var.units = 'see source code' + else: + Var.units = 'fraction' + +# global attributes + version = '1.0' + ncFile.Author = 'TUV Data Xformer ' + version + now = dt.now() + ncFile.created = now.strftime("%Y-%m-%d %H:%M:%S") + if( var_typ == 'cross_section' ): + ncFile.title = 'Cross section parameters' + else: + ncFile.title = 'Quantum yield parameters' + ncFile.file = InpFileSpecs + +""" +Transform ascii data file(s) to netCDF counterpart +""" +def xform_to_netCDF(nFile,phtDictionary,ncd_path): + +# cross section + if( 'cross-sections' in phtDictionary ): +# form netCDF file for cross sections + molecule = phtDictionary['molecule'] + + print(f'\nProcessing {molecule} cross section') + xsects = phtDictionary['cross-sections'] + nxsects = len(xsects) + print(f' There are {nxsects} cross section files') + dataTray = [] + interpolationTemps = [] + Headers = [] + +# loop over ascii input data files + for xsect in xsects: + ncdFilespec = ncd_path + '/' + molecule + '_cross_section_' + str(nFile) + '.nc' +# create the netcdf dataset + print(f'\nCreating netCDF file {ncdFilespec}') + + try: + ncFile = Dataset(ncdFilespec,mode='w',format='NETCDF4_CLASSIC') + except: + print(f'Failed to create netCDF4 dataset {ncdFilespec}') + sys.exit(-1) + +# 1st data column wavelength grid? + if( 'has lambda grid' in xsect ): + hasLambdaGrid = xsect['has lambda grid'] + else: + hasLambdaGrid = True + +# interpolation temperatures? + if( 'interpolation temperature' in xsect ): + interpolationTemps = xsect['interpolation temperature'] + print("\nInterpolation temperatures:") + print(interpolationTemps) + + header = '' + header = read_data_file(xsect,dataTray) + Headers.append(header) + + InpFileSpecs = xsect['filespec'] + + print(f'\nThere are {len(dataTray)} arrays in dataTray') + print(f'\nShape data array is {dataTray[0].shape}') + stuff_netCDF_file(ncFile,interpolationTemps,dataTray,hasLambdaGrid,InpFileSpecs,'cross_section',Headers) +# close netcdf file + ncFile.close() + + print('\n') + +# quantum yield + if( 'quantum-yields' in phtDictionary ): +# form netCDF file for cross sections + molecule = phtDictionary['molecule'] + + print(f'\nProcessing {molecule} quantum yield') + qylds = phtDictionary['quantum-yields'] + nqylds = len(qylds) + print(f' There are {nqylds} quantum yield files') + dataTray = [] + interpolationTemps = [] + Headers = [] + +# loop over ascii input data files + for qyld in qylds: + ncdFilespec = ncd_path + '/' + molecule + '_quantum_yield_' + str(nFile) + '.nc' +# create the netcdf dataset + print(f'\nCreating netCDF file {ncdFilespec}') + + try: + ncFile = Dataset(ncdFilespec,mode='w',format='NETCDF4_CLASSIC') + except: + print(f'Failed to create netCDF4 dataset {ncdFilespec}') + sys.exit(-1) + +# 1st data column wavelength grid? + if( 'has lambda grid' in qyld ): + hasLambdaGrid = qyld['has lambda grid'] + else: + hasLambdaGrid = True + +# interpolation temperatures? + if( 'interpolation temperature' in qyld ): + interpolationTemps = qyld['interpolation temperature'] + print("\nInterpolation temperatures:") + print(interpolationTemps) + + header = '' + header = read_data_file(qyld,dataTray) + Headers.append(header) + + InpFileSpecs = qyld['filespec'] + + print(f'\nThere are {len(dataTray)} arrays in dataTray') + print(f'\nShape data array is {dataTray[0].shape}') + stuff_netCDF_file(ncFile,interpolationTemps,dataTray,hasLambdaGrid,InpFileSpecs,'quantum_yield',Headers) +# close netcdf file + ncFile.close() From cb6efa67f135f86b508e11fd0fecc9e9426abd0d Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Thu, 11 Jan 2024 11:14:19 -0800 Subject: [PATCH 2/2] fix memory leak in tests --- test/unit/tuv_doug/data_sets.F90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/unit/tuv_doug/data_sets.F90 b/test/unit/tuv_doug/data_sets.F90 index 50874616..f62da19c 100644 --- a/test/unit/tuv_doug/data_sets.F90 +++ b/test/unit/tuv_doug/data_sets.F90 @@ -89,7 +89,7 @@ subroutine test_data_set( ) default = 1.0e-6_dk ) call config_pair%get( "mask", mask_points_config, Iam, & found = found ) - if (found) then + if( found ) then mask_points_iter => mask_points_config%get_iterator( ) allocate( mask_points( mask_points_config%number_of_children( ) ) ) do i = 1, size( mask_points ) @@ -99,8 +99,9 @@ subroutine test_data_set( ) call mask_point_config%get( "index", mask_points( i ), Iam ) end do call assert( 888375064, .not. mask_points_iter%next( ) ) + deallocate( mask_points_iter ) else - allocate(mask_points(0)) + allocate( mask_points(0) ) end if ! Load and test cross section