[ncl-talk] Reversing direction of integration yields different values

Thu Apr 5 11:28:04 MDT 2018

Hi Andreas,

Did you get this figured out? Alan provided some good suggestions about
avoiding do loops, so hopefully you were able to incorporate those.

The other issue is that you were looping across "k", but starting at index
1:

; #1 Method
do j = 0, nlat-1
  do k = 1, nlev-1
    psistar(0,:) = 0
    if (k .eq. 0) then
      psistar(k,j) = 2*pi*r*COSLAT(j)*(VSTAR(k,j)*dp(k))/g
    else if (k .gt. 0)
      psistar(k,j) = psistar(k-1,j) +
(2*pi*r*COSLAT(j)*(VSTAR(k,j)*dp(k)))/g
    end if
  end if
end do
end do

You have a check in the loop for k.eq.0, but this will never get reached
because k starts at 1.

If you are still having problems with the script, please post back to
ncl-talk and include the latest version of your code.

Thanks,

--Mary

On Tue, Mar 27, 2018 at 2:38 AM, Andreas Chrysanthou <eeac at leeds.ac.uk>
wrote:

> Hi Alan,
>
> Thanks for this, both for the suggestions and the effort to debug the code.
>
> 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.
>
> So unfortunately I’m back to square 1.
>
> Cheers,
> Andreas
>
>
> On 27 Mar 2018, at 04:56, Alan Brammer <abrammer at albany.edu> wrote:
>
>
> I think the problem may be here:
> if (k .eq. 0) then
>     dp(k) = lev(k)*100
>
> So in one direction dp(0)  = 0.1 *100
> in the other direction dp(0) = 1000*100
>
> That seems like it would introduce some large differences depending on the
> direction of the array.
>
> There is this function to calculate pressure differences
> https://www.ncl.ucar.edu/Document/Functions/Built-in/dpres_plevel.shtml
>
>
> Also, some minor suggestions:
> 1.
> I think this can be broken out of a loop and just computed as an array.
> Should be marginally faster in NCL.
> COSLAT(j) = cos(lat(j)*d2rad)
> COSLAT = cos(lat*d2rad)
>
>
> 2.
> in Method 1.
> 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.
>
> do k=0...
> ....
> psistar(k,:) = 2*pi*r*COSLAT*(VSTAR(k,:)*dp(k))/g
> ....
> end do
>
>
>
> On Mon, Mar 26, 2018 at 11:19 AM, Andreas Chrysanthou <eeac at leeds.ac.uk>
> wrote:
>
>> Hi NCL users,
>>
>> I came up into something really weird unless I’m ignoring something
>> really big, which is not likely.
>>
>> 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).
>>
>> 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?
>>
>> Attached is a snip of my script with the relevant part. Also two plots of
>> both version to compare the differences.
>>
>> ; Constants
>> pi    = 4.*atan(1.)
>> r     = 6371*10^3   ; Earth Radius
>> g     = 9.81         ; m/s2 Gravity
>> d2rad = 0.017453  ; degrees to radians
>>
>> varname = "vstar"
>>
>> vstar = ff1->$varname$(0:599,:,:,0) ; normal array pressure levels
>> (1000->01.hPa)
>> ;vstar = vstar(:,::-1,:) ; uncomment if integrating from top -> bottom
>>
>> VSTAR = dim_avg_n_Wrap(vstar, 0) ; just a timemean
>> printVarSummary(VSTAR)
>> printMinMax(VSTAR,1)
>>
>> lat = VSTAR&lat
>> lev = VSTAR&lev
>>
>> nlev = dimsizes(lev)
>> nlat = dimsizes(lat)
>>
>> dp = new(nlev,typeof(lev))
>> do k = 0,nlev-1
>>   if (k .eq. 0) then
>>     dp(k) = lev(k)*100
>>   else
>>     dp(k) = abs((lev(k)-lev(k-1))*100)
>>   end if
>> end do
>>
>> COSLAT = new(nlat,typeof(lat))
>> do j = 0, nlat-1
>>   COSLAT(j) = cos(lat(j)*d2rad)
>> end do
>> print(COSLAT)
>>
>> psistar   = new((/nlev,nlat/),"double")      ; allocate space
>> psistar!0 = "lev"
>> psistar!1 = "lat"
>> psistar&lev = lev
>> psistar&lat = lat
>>
>> ;***********************************************************
>> *********************************************
>> ; NOTE these two methods yield identical results, correct if lev axis is
>> reversed. When not, integrating
>> ; bottom to the top yields unrealistic results with both formulations
>>
>> ; #1 Method
>> do j = 0, nlat-1
>>   do k = 1, nlev-1
>>     psistar(0,:) = 0
>>     if (k .eq. 0) then
>>       psistar(k,j) = 2*pi*r*COSLAT(j)*(VSTAR(k,j)*dp(k))/g
>>     else if (k .gt. 0)
>>       psistar(k,j) = psistar(k-1,j) + (2*pi*r*COSLAT(j)*(VSTAR(k,j)*
>> dp(k)))/g
>>     end if
>>   end if
>> end do
>> end do
>>
>> ; #2 Method
>> ;do j = 0, nlat-1
>> ;  do k = 0 , nlev-1
>> ;    ;psistar(0,j) = 0
>> ;    psistar(k,j) = 2*pi*r*COSLAT(j)*dim_sum(VSTAR(0:k,j)*dp(0:k))/g
>> ;    psistar(0,j) = 0
>> ;    ;psistar(nlev-1,j) = 0
>> ;  end do
>> ;end do
>>
>> ;***********************************************************
>> ************************************************
>>
>>
>>
>> Cheers,
>> Andreas
>>
>>
>>
>> _______________________________________________
>> ncl-talk mailing list
>> ncl-talk at ucar.edu
>> List instructions, subscriber options, unsubscribe:
>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
>>
>
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180405/35010a71/attachment.html>