[ncl-talk] grad_latlon_cfd problem

[Huan Wang]

Hi,
In regard to this question,

I'm trying to use grad_latlon_cfd function by loading contributed.ncl for version 6.
My codes are here,

;NCL tutorial script: vert_1.ncl
;************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/skewt_func.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "/home/wangh1/contributed.ncl"
;************************************************

begin

;************************************************
; file handling
;************************************************
fn  =systemfunc("ls /*/con.cam.h0.*.nc"); define filename, 20 years summer files
in  = addfiles(fn,"r")                         ; open netcdf file
ListSetType(in, "join")

;************************************************
; read needed variables from file
;************************************************
PS   = in[:]->PS                                ; get pressure
lat= in[1]->lat
lon=in[1]->lon
gradLatLon= grad_latlon_cfd(PS,PS&lat,PS&lon,True,False)
glon= gradLatLon[0]
glat= gradLatLon[1]
end


And the result is:
fatal:Assignment type mismatch, right hand side can't be coerced to type of left hand side
fatal:["Execute.c":8578]:Execute: Error occurred at or near line 58 in file /home/wangh1/contributed.ncl

fatal:["Execute.c":8578]:Execute: Error occurred at or near line 32 in file bou2.ncl

Contributed.ncl file I used is:

undef("grad_latlon_cfd")
function grad_latlon_cfd(z:numeric, lat[*]:numeric, lon[*]:numeric, rCyclic[1]:logical, opt[1]:logical)
local dimz, rankz, nlat, mlon, rad, re, con, dlon, dx, cgx, cgy, gradLatLon
begin
  dimz  = dimsizes(z)
  rankz = dimsizes(dimz)

  if (rankz.lt.2 .or. rankz.gt.4) then
      print("grad_latlon_cfd: illegal rank: rankz="+rankz)
      exit
  end if

  mlon  = dimz(rankz-1)  
  nlat  = dimz(rankz-2) 
;************************************************
; Miscellaneous
;************************************************
  rad   = 4.0*atan(1.0)/180.0
  re    = 6.37122e6  
  con   = re*rad                 ; one deg lat = 111198.8 meters

;************************************************
; Use 'center_finite_diff_n for meridional (Y) gradient
;************************************************
  cgy  = center_finite_diff_n( z, lat, False ,0,rankz-2) 
  cgy  = cgy/con
  copy_VarCoords( z, cgy)                         ; add mets 
  cgy at long_name = "cfd: meridional gradient"
  cgy at units     = "?/m"

;************************************************
; Use cfd for zonal (X) gradients
; Generally, these are much smaller than the meridional (Y) gradients
; This assumes that the longitudes are equally spaced
; Pre-allocate space for gradients
;************************************************

  dlon = (lon(2)-lon(1))*0.0174533   ; convert to radians
  dx   = re*cos(0.0174533*lat)*dlon  ; different at each latitude; dx(nlat)
  cgx  = new( dimz, typeof(z), getFillValue(z) )  ; lon=>X
;;       These don't work
;;DX   = conform_dims( dimz, dx, rankz-2)
;;CGX  = center_finite_diff_n (z, DX, rCyclic,0,0)
;;copy_VarCoords( z, CGX)            
;;CGX at long_name = "gradx: ZONAL GRADIENT"
;;CGX at units     = "?/m"

  if (rankz.eq.2) then
      do nl=0,nlat-1                     ; loop over each latitude
          cgx(nl,:)     = center_finite_diff_n (z(nl,:), dx(nl), rCyclic,0,0)
      end do
  else if (rankz.eq.3) then
      do nl=0,nlat-1                     ; loop over each latitude
          cgx(:,nl,:)   = center_finite_diff_n (z(:,nl,:), dx(nl), rCyclic,0,0)
      end do
  else if (rankz.eq.4) then
      do nl=0,nlat-1                     ; loop over each latitude
          cgx(:,:,nl,:) = center_finite_diff_n (z(:,:,nl,:), dx(nl), rCyclic,0,0)
      end do
    end if      ; rankz=4
   end if       ; rankz=3
  end if        ; rankz=2

  copy_VarCoords( z, cgx)            
  cgx at long_name = "cfd: zonal gradient"
  cgx at units     = "?/m"

  gradLatLon = [/ cgy, cgx /]   ; return two variables as a type 'list'
  return( gradLatLon )
end



Thank you!


Huan