;************************************************
;This script calculates the temperature gradient
;************************************************

begin
;************************************************
; Variable file handling
;************************************************
  f  = addfile("thetae_ncep_jan6h.nc","r")       ; open netcdf file
  T   = f->theta_e(62,{850},:,:)


  T = T(::-1,:)
  T = lonPivot(T, 30)

;************************************************
; Miscellaneous
;************************************************

  dimT = dimsizes(T)
  nlat  = dimT(0)                       ; 94
  mlon  = dimT(1)                       ; 192

  rad   = 4.*atan(1.)/180.
  lat   = T&lat
  lon   = T&lon

  re    = 6.37122e6                      ; spherical earth
  con   = re*rad                         ; one deg lat = 111198.8 meters
  scly  = 1e5                            ; scale => nicer plots
  sclx  = 1e5

;************************************************
; SPHERICAL HARMONICS:
; Use "highly accurate" spherical harmonics procedure (gradsg)
; to compute zonal (X) and meridional (Y) gradients.
;************************************************
                                         ; pre-allocate space for return gradients
  TGX_gradsg = new( dimT, typeof(T), getFillValue(T) )  ; lon=>X
  TGY_gradsg = new( dimT, typeof(T), getFillValue(T) )  ; lat=>Y
  gradsg(T, TGX_gradsg, TGY_gradsg)     ; procedure for gaussian grids

  copy_VarCoords(T, TGX_gradsg)         ; add meta data
  copy_VarCoords(T, TGY_gradsg) 
  TGX_gradsg@long_name = "TGX: gradsg"
  TGX_gradsg@units     = "K/s"
  TGY_gradsg@long_name = "TGY: gradsg"
  TGY_gradsg@units     = "K/s"

;**********************************************
;Plot
;**********************************************
  wks = gsn_open_wks("png","gradient")          ; send graphics to PNG file

  res              = True
  res@gsnMaximize  = True     ; maximize plot in frame
  res@cnFillOn     = True     ; turn on contour fill
  res@cnLinesOn    = False    ; turn off contour fill
  res@tiMainString = "Test"

  res@mpMinLatF = -5
  res@mpMaxLatF = 50
  res@mpMinLonF = 100
  res@mpMaxLonF = 200

  plot = gsn_csm_contour_map(wks,TGY_gradsg,res)
end