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

[Andreas Chrysanthou]

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<mailto: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<mailto: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<mailto: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/20180327/ae39a0ce/attachment.html>