<div dir="ltr"><div class="gmail_default" style="font-family:verdana,sans-serif"><br></div><div class="gmail_default" style="font-family:verdana,sans-serif">I think the problem may be here:</div><div class="gmail_default" style="font-family:verdana,sans-serif"><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px">if (k .eq. 0) then<br></div><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px"> dp(k) = lev(k)*100</div><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px"><br></div><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px">So in one direction dp(0) = 0.1 *100</div><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px">in the other direction dp(0) = 1000*100</div><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px"><br></div><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px">That seems like it would introduce some large differences depending on the direction of the array. </div><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px"><br></div><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px">There is this function to calculate pressure differences <a href="https://www.ncl.ucar.edu/Document/Functions/Built-in/dpres_plevel.shtml">https://www.ncl.ucar.edu/Document/Functions/Built-in/dpres_plevel.shtml</a></div><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px"><br></div><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px"><br></div><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px">Also, some minor suggestions:</div><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px">1.</div><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px">I think this can be broken out of a loop and just computed as an array. Should be marginally faster in NCL. </div><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px">COSLAT(j) = cos(lat(j)*d2rad)<br></div><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px">COSLAT = cos(lat*d2rad)<br></div><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px"><br></div><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px"><br></div><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px">2. </div><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px">in Method 1. </div><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px">The j, loop can be removed and just compute on all lats at once. Nested do loops become slow rather quickly in NCL, better to vectorise as mush as possible.</div><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px"><br></div><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px">do k=0...</div><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px">....</div><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px">psistar(k,:) = 2*pi*r*COSLAT*(VSTAR(k,:)*<wbr>dp(k))/g<br></div><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px">....</div><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px">end do</div><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px"><br></div><div style="color:rgb(0,0,0);font-family:arial,sans-serif;font-size:12.800000190734863px"><br></div></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Mon, Mar 26, 2018 at 11:19 AM, Andreas Chrysanthou <span dir="ltr"><<a href="mailto:eeac@leeds.ac.uk" target="_blank">eeac@leeds.ac.uk</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div>
<div style="word-wrap:break-word">Hi NCL users,
<div><br>
</div>
<div>I came up into something really weird unless I’m ignoring something really big, which is not likely. </div>
<div><br>
</div>
<div>I’m trying to calculate the mass streamfunction (psistar) from the meridional velocity and when I’m integrating downwards (reversing the v wind array) from Top of the Atmosphere (TOA - 0.1hPa) to Bottom (1000hPa) the calculation yields some very
realistic values and the proper structure. (Magnitude 10^9 Kg s-1). </div>
<div><br>
</div>
<div>If I try to do that the other way round (from bottom-TOA) I’m getting a difference of 10^3 in the magnitude and not a very smooth structure and this is really counter intuitive. Any ideas If I’m doing something wrong?</div>
<div><br>
</div>
<div>Attached is a snip of my script with the relevant part. Also two plots of both version to compare the differences.</div>
<div><br>
</div>
<div>
<div>; Constants </div>
<div>pi = 4.*atan(1.)</div>
<div>r = 6371*10^3 ; Earth Radius</div>
<div>g = 9.81 ; m/s2 Gravity</div>
<div>d2rad = 0.017453 ; degrees to radians</div>
<div><br>
</div>
<div>varname = "vstar"</div>
<div><br>
</div>
<div>vstar = ff1->$varname$(0:599,:,:,0) ; normal array pressure levels (1000->01.hPa)</div>
<div>;vstar = vstar(:,::-1,:) ; uncomment if integrating from top -> bottom</div>
<div><br>
</div>
<div>VSTAR = dim_avg_n_Wrap(vstar, 0) ; just a timemean</div>
<div>printVarSummary(VSTAR)</div>
<div>printMinMax(VSTAR,1)</div>
<div><br>
</div>
<div>lat = VSTAR&lat</div>
<div>lev = VSTAR&lev</div>
<div><br>
</div>
<div>nlev = dimsizes(lev)</div>
<div>nlat = dimsizes(lat)</div>
<div><br>
</div>
<div>dp = new(nlev,typeof(lev))</div>
<div>do k = 0,nlev-1</div>
<div> if (k .eq. 0) then</div>
<div> dp(k) = lev(k)*100</div>
<div> else</div>
<div> dp(k) = abs((lev(k)-lev(k-1))*100)</div>
<div> end if</div>
<div>end do</div>
<div><br>
</div>
<div>COSLAT = new(nlat,typeof(lat))</div>
<div>do j = 0, nlat-1</div>
<div> COSLAT(j) = cos(lat(j)*d2rad)</div>
<div>end do</div>
<div>print(COSLAT)</div>
<div><br>
</div>
<div>psistar = new((/nlev,nlat/),"double") ; allocate space</div>
<div>psistar!0 = "lev"</div>
<div>psistar!1 = "lat"</div>
<div>psistar&lev = lev</div>
<div>psistar&lat = lat</div>
<div><br>
</div>
<div>;*****************************<wbr>******************************<wbr>******************************<wbr>***************</div>
<div>; NOTE these two methods yield identical results, correct if lev axis is reversed. When not, integrating</div>
<div>; bottom to the top yields unrealistic results with both formulations</div>
<div><br>
</div>
<div>; #1 Method</div>
<div>do j = 0, nlat-1</div>
<div> do k = 1, nlev-1 </div>
<div> psistar(0,:) = 0</div>
<div> if (k .eq. 0) then</div>
<div> psistar(k,j) = 2*pi*r*COSLAT(j)*(VSTAR(k,j)*<wbr>dp(k))/g</div>
<div> else if (k .gt. 0)</div>
<div> psistar(k,j) = psistar(k-1,j) + (2*pi*r*COSLAT(j)*(VSTAR(k,j)*<wbr>dp(k)))/g</div>
<div> end if</div>
<div> end if</div>
<div>end do</div>
<div>end do</div>
<div><br>
</div>
<div>; #2 Method</div>
<div>;do j = 0, nlat-1</div>
<div>; do k = 0 , nlev-1</div>
<div>; ;psistar(0,j) = 0</div>
<div>; psistar(k,j) = 2*pi*r*COSLAT(j)*dim_sum(<wbr>VSTAR(0:k,j)*dp(0:k))/g</div>
<div>; psistar(0,j) = 0</div>
<div>; ;psistar(nlev-1,j) = 0 </div>
<div>; end do</div>
<div>;end do</div>
<div><br>
</div>
<div>;*****************************<wbr>******************************<wbr>******************************<wbr>******************</div>
<div><br>
</div>
<div><br>
</div>
<div></div>
</div>
</div>
<div style="word-wrap:break-word">
<div>
<div></div>
</div>
</div>
<div style="word-wrap:break-word">
<div>
<div></div>
<div><br>
</div>
<div>Cheers,</div>
<div>Andreas</div>
<div>
<div style="color:rgb(0,0,0);letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;word-wrap:break-word">
<div style="color:rgb(0,0,0);letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;word-wrap:break-word">
<div style="color:rgb(0,0,0);letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;word-wrap:break-word">
<div style="color:rgb(0,0,0);letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;word-wrap:break-word">
<div style="color:rgb(0,0,0);letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;word-wrap:break-word">
<div><br>
</div>
<div></div>
</div>
</div>
</div>
</div>
</div>
</div>
<br>
</div>
</div>
</div>
<br>______________________________<wbr>_________________<br>
ncl-talk mailing list<br>
<a href="mailto:ncl-talk@ucar.edu">ncl-talk@ucar.edu</a><br>
List instructions, subscriber options, unsubscribe:<br>
<a href="http://mailman.ucar.edu/mailman/listinfo/ncl-talk" rel="noreferrer" target="_blank">http://mailman.ucar.edu/<wbr>mailman/listinfo/ncl-talk</a><br>
<br></blockquote></div><br></div>