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

[Andreas Chrysanthou]

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


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180326/81ec8c93/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: bottom_top.pdf
Type: application/pdf
Size: 41621 bytes
Desc: bottom_top.pdf
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180326/81ec8c93/attachment.pdf>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: top_bottom.pdf
Type: application/pdf
Size: 51415 bytes
Desc: top_bottom.pdf
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180326/81ec8c93/attachment-0001.pdf>