The result of your pow() function has no solution in the reals (the result is complex). To get the behavior you want do:<div><br /></div><div>pow(a,k) = a*pow(abs(a),k-1)</div><div><br /></div><div>or something similar.</div><div><br /></div><div>-Dave<br /><br /><span>On 01/18/13, <b class="name">Jeff Duda </b> &lt;jeffduda319@gmail.com&gt; wrote:</span><blockquote cite="mid:CAAig09DntUuetXdOGho8Fqk3OCYZY8HpVEX17yu26RA=9xyc=A@mail.gmail.com" class="iwcQuote" style="border-left: 1px solid #00F; padding-left: 13px; margin-left: 0;" type="cite"><div class="mimepart text html">Hello,<br />I am trying to compute some thermodynamic parameters, shown below.<br clear="all" /><br /> 'define e = rhprs*6.11*exp((2.5e6/461.5)*((1/273)-(1/tmpprs)))'<br /> 'define r = 622*e/(lev-e)'<br /> 'define condtemp = 55 + 1 / (1/(tmpprs-55) - log(rhprs/100)/2840)'<br />
 'define thetaDL = tmpprs*pow(1000/(lev-e),0.2854)*pow(tmpprs/condtemp,0.0028*r)'<br /> 'define thetae = thetaDL*exp(((3.036/condtemp)-0.00178)*r*(1+0.448+r/1000))'<br /> 'define thetaw = 45.114 - 51.489*pow(thetae/273.15,-3.504)'<br />
 'define LSI = max(thetaw,lev=1000,lev=500) - max(thetaw,lev=1000,lev=900)'<br /><br />I get a lot of NaNs.  They are coming from this line:<br /> 'define thetaDL = tmpprs*pow(1000/(lev-e),0.2854)*pow(tmpprs/condtemp,0.0028*r)'<br />
<br />And in fact, from this piece:<br />pow(1000/(lev-e),0.2854)<br /><br />I know this because I tested this for a model grid column (from 1000 - 500 mb). <br />ga-&gt; d 1000/(lev-e)            <br />Printing Grid -- 13 Values -- Undef = -9.99e+08<br />
-3.55595 -4.04332 -4.69364 -9.98891 39.3747 -5.3341 -11.5849 4.19346 <br />2.07427 2.42358 2.11214 2.1287 2.24423 <br />ga-&gt; d pow(1000/(lev-e),0.2854)<br />Printing Grid -- 13 Values -- Undef = -9.99e+08<br />-nan -nan -nan -nan 2.85282 -nan -nan 1.5055 <br />
1.2315 1.28743 1.23787 1.24063 1.25949 <br /><br />For whatever reason, the pow() function doesn't like negative bases, because for example, I can use a calculator to compute -3.55595^0.2854 as -1.43629.  Why, then, does the pow() function give a -nan as a result?<br />
<br />Jeff Duda<br /><br />-- <br />Jeff Duda<br />Graduate research assistant<br />University of Oklahoma School of Meteorology<br />Center for Analysis and Prediction of Storms<br />
</div></blockquote></div>