<div dir="ltr"><div class="gmail_default" style="font-size:small">Hi Andreas,</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">Did you get this figured out? Alan provided some good suggestions about avoiding do loops, so hopefully you were able to incorporate those.</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">The other issue is that you were looping across "k", but starting at index 1:</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default"><div class="gmail_default"><div class="gmail_default">; #1 Method</div><div class="gmail_default">do j = 0, nlat-1</div><div class="gmail_default">  do k = 1, nlev-1 </div><div class="gmail_default">    psistar(0,:) = 0</div><div class="gmail_default">    if (k .eq. 0) then</div><div class="gmail_default">      psistar(k,j) = 2*pi*r*COSLAT(j)*(VSTAR(k,j)*dp(k))/g</div><div class="gmail_default">    else if (k .gt. 0)</div><div class="gmail_default">      psistar(k,j) = psistar(k-1,j) + (2*pi*r*COSLAT(j)*(VSTAR(k,j)*dp(k)))/g</div><div class="gmail_default">    end if</div><div class="gmail_default">  end if</div><div class="gmail_default">end do</div><div class="gmail_default">end do</div><div class="gmail_default"><br></div><div class="gmail_default">You have a check in the loop for k.eq.0, but this will never get reached because k starts at 1.</div><div class="gmail_default"><br></div><div class="gmail_default">If you are still having problems with the script, please post back to ncl-talk and include the latest version of your code.</div><div class="gmail_default"><br></div><div class="gmail_default">Thanks,</div><div class="gmail_default"><br></div><div class="gmail_default">--Mary</div><div class="gmail_default"><br></div></div></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Mar 27, 2018 at 2:38 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 style="word-wrap:break-word">
Hi Alan,
<div><br>
</div>
<div>Thanks for this, both for the suggestions and the effort to debug the code.</div>
<div><br>
</div>
<div>The first element of the dp unfortunately is not playing any part to it, as I’ve tried it with either dp(0) =-dp(1) or with just setting it as zero. </div>
<div><br>
</div>
<div>So unfortunately I’m back to square 1.</div>
<div><br>
</div>
<div>Cheers,</div>
<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>Andreas</div>
<div><br>
</div>
<div></div>
</div>
</div>
</div>
</div>
</div>
</div><div><div class="h5">
<br>
<div>
<blockquote type="cite">
<div>On 27 Mar 2018, at 04:56, Alan Brammer <<a href="mailto:abrammer@albany.edu" target="_blank">abrammer@albany.edu</a>> wrote:</div>
<br class="m_3952703453369567090Apple-interchange-newline">
<div>
<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="font-family:arial,sans-serif;font-size:12.800000190734863px">
if (k .eq. 0) then<br>
</div>
<div style="font-family:arial,sans-serif;font-size:12.800000190734863px">
    dp(k) = lev(k)*100</div>
<div style="font-family:arial,sans-serif;font-size:12.800000190734863px">
<br>
</div>
<div style="font-family:arial,sans-serif;font-size:12.800000190734863px">
So in one direction dp(0)  = 0.1 *100</div>
<div style="font-family:arial,sans-serif;font-size:12.800000190734863px">
in the other direction dp(0) = 1000*100</div>
<div style="font-family:arial,sans-serif;font-size:12.800000190734863px">
<br>
</div>
<div style="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="font-family:arial,sans-serif;font-size:12.800000190734863px">
<br>
</div>
<div style="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" target="_blank">https://www.ncl.ucar.edu/<wbr>Document/Functions/Built-in/<wbr>dpres_plevel.shtml</a></div>
<div style="font-family:arial,sans-serif;font-size:12.800000190734863px">
<br>
</div>
<div style="font-family:arial,sans-serif;font-size:12.800000190734863px">
<br>
</div>
<div style="font-family:arial,sans-serif;font-size:12.800000190734863px">
Also, some minor suggestions:</div>
<div style="font-family:arial,sans-serif;font-size:12.800000190734863px">
1.</div>
<div style="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="font-family:arial,sans-serif;font-size:12.800000190734863px">
COSLAT(j) = cos(lat(j)*d2rad)<br>
</div>
<div style="font-family:arial,sans-serif;font-size:12.800000190734863px">
COSLAT = cos(lat*d2rad)<br>
</div>
<div style="font-family:arial,sans-serif;font-size:12.800000190734863px">
<br>
</div>
<div style="font-family:arial,sans-serif;font-size:12.800000190734863px">
<br>
</div>
<div style="font-family:arial,sans-serif;font-size:12.800000190734863px">
2. </div>
<div style="font-family:arial,sans-serif;font-size:12.800000190734863px">
in Method 1. </div>
<div style="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="font-family:arial,sans-serif;font-size:12.800000190734863px">
<br>
</div>
<div style="font-family:arial,sans-serif;font-size:12.800000190734863px">
do k=0...</div>
<div style="font-family:arial,sans-serif;font-size:12.800000190734863px">
....</div>
<div style="font-family:arial,sans-serif;font-size:12.800000190734863px">
psistar(k,:) = 2*pi*r*COSLAT*(VSTAR(k,:)*dp(<wbr>k))/g<br>
</div>
<div style="font-family:arial,sans-serif;font-size:12.800000190734863px">
....</div>
<div style="font-family:arial,sans-serif;font-size:12.800000190734863px">
end do</div>
<div style="font-family:arial,sans-serif;font-size:12.800000190734863px">
<br>
</div>
<div style="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)*d<wbr>p(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(VSTAR<wbr>(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="letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;word-wrap:break-word">
<div style="letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;word-wrap:break-word">
<div style="letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;word-wrap:break-word">
<div style="letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;word-wrap:break-word">
<div style="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" target="_blank">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/mailma<wbr>n/listinfo/ncl-talk</a><br>
<br>
</blockquote>
</div>
<br>
</div>
</div>
</blockquote>
</div>
<br>
</div></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>