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

[Alan Brammer]

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
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180326/c5444e50/attachment.html>