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

Fri Apr 13 12:00:17 MDT 2018

Hi Mary,

Sorry for the late reply, I was offline.

Unfortunately I didn’t solve it yet, the k index was wrongly set to 1 when I last sent the script, it was set at 0 initially. So I still can’t understand why I’m getting two orders of magnitude difference in the results when I integrate from the surface to the top. Attached is the latest version of my script, not that different from the previous one.

Any kind of help or pointers will be much appreciated!

vstar = ff1->$varname$(0:599,:,:,0)
;vstar = vstar(:,::-1,:) uncomment if integrating from top -> bottom

***************************************************
***************************************************
VSTAR = dim_avg_n_Wrap(vstar, 0)
printVarSummary(VSTAR)
printMinMax(VSTAR,1)
delete(vstar)

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)

;******************************************
; Simplified Version of the above
;COSLAT = cos(lat*d2rad)
;******************************************
psistar   = new((/nlev,nlat/),"double")      ; allocate space
psistar!0 = "lev"
psistar!1 = "lat"
psistar&lev = lev
psistar&lat = lat
printVarSummary(psistar)
print(lev)

;********************************************************************************************************
; NOTE these two methods yield identical results, correct if lev axis is reversed. When not, integrating
; bottom to the top yields unrealistic results

 #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 (Use of built in function of NCL - dim_sum)
;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_v(0:k))/g
;    psistar(0,j) = 0
;    ;psistar(nlev-1,j) = 0
;  end do
;end do

Cheers,

Andreas

On 5 Apr 2018, at 19:28, Mary Haley <haley at ucar.edu<mailto:haley at ucar.edu>> wrote:

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<mailto: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<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

_______________________________________________
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/20180413/3d86e2f7/attachment.html>