;   Script to produce plots for a WRF real-data run,
;   with the ARW coordinate dynamics option.

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "./WRFUserARW_add.ncl"

begin

; We generate plots, but what kind do we prefer?
  type = "x11"
; type = "pdf"
; type = "ps"
; type = "ncgm"
  wks = gsn_open_wks(type,"plt_addfiles")

; Set some basic resources
  res = True
  res@MainTitle                   = "REAL-TIME WRF"

  pltres = True
  mpres = True

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
; Make a list of all files we are interested in
  DATADir = "./"
  ;FILES = systemfunc ("csh -c ' ls -1 " + DATADir + "wrfout?.nc'")
  FILES = systemfunc ("csh -c ' ls -1 " + DATADir + "wrfout_new?.nc'")
  ;FILES = systemfunc ("csh -c ' ls -1 " + DATADir + "wrfout_d01_2000-*00'")

  numFILES = dimsizes(FILES)
  ;do ff = 0,numFILES-1
  ;  FILES(ff) = FILES(ff)+".nc"
  ;end do

  print("numFILES = " + numFILES)
  print(FILES)

  f = addfiles(FILES,"r")
  a = addfile("./wrfout_d01_2000-01-24_12:00:00.nc","r")

  var1 = wrf_user_getvar(a,"T",3)                                 ; time 3 from a file with 5 times
  var2 = wrf_user_getvar_from_files(FILES,"T",-1,True)            ; all times from split files (total 5 times)
  var3 = wrf_user_getvar_from_files(FILES,"T",(/1,4,2/),True)     ; get times 1,3 from split files (from 5 times avail)
  var4 = wrf_user_getvar_from_files(FILES,"T",3,True)             ; time 3 from split files (total of 5 time avail)
  times = wrf_user_getvar_from_files(FILES,"times",(/0,4,1/),True) ; get all times
  print(times)


  var1@description = "Theta for var1"
  var2@description = "Theta for var2"
  var3@description = "Theta for var3"
  var4@description = "Theta for var4"

  ; Plotting options
  res@TimeLabel = times(3)    ; I happend to know that this match all the plots below
  opts = res
  opts@cnFillOn = True

  contour_var1 = wrf_contour(a,wks,var1(0,:,:),opts)             ; only have one time - plot it
  contour_var2 = wrf_contour(f[0],wks,var2(3,0,:,:),opts)        ; have all the times, plot 3rd
  contour_var3 = wrf_contour(f[0],wks,var3(1,0,:,:),opts)        ; only got 1,3, so we want to plot 1
  contour_var4 = wrf_contour(f[0],wks,var4(0,0,:,:),opts)        ; only got the 3rd - so plot it

  plot = wrf_map_overlays(a,wks,(/contour_var1/),pltres,mpres)
  plot = wrf_map_overlays(f[0],wks,(/contour_var2/),pltres,mpres)
  plot = wrf_map_overlays(f[0],wks,(/contour_var3/),pltres,mpres)
  plot = wrf_map_overlays(f[0],wks,(/contour_var4/),pltres,mpres)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

end