;   Example 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 "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"   

external CFTD "./cftd.so"

begin
;
; The WRF ARW input file.  
; This needs to have a ".nc" appended, so just do it.

  filesa = systemfunc ("ls /glade/scratch/yingsong/WRFV3.6.1_real_out_new/Aug27/wrfout_00Z-21Z_control_r0.08_15min/wrfout_d02_2009-08-27*") ; file paths
    a = addfiles(filesa,"r")
    ListSetType (a, "cat")

; type = "pdf"
 type = "ps"
  wks = gsn_open_wks(type,"plt_wa_time_evo_Aug27_test")
  gsn_define_colormap(wks,"tbrAvg1")      ; choose colormap


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

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Which times and how many time steps are in the data set?
  FirstTime = True
  times = wrf_user_getvar(a,"times",-1)  ; get all times in the file
  ntimes = dimsizes(times)         ; number of times in the file

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   wa =wrf_user_getvar(a, "wa", -1)
   wa@_FillValue = -9999.0
   z =wrf_user_getvar(a, "z",-1)

max_wa = max(wa)
print("max_wa ="+max_wa)
min_wa = min(wa)
print("min_wa ="+min_wa)

  dsizes_x = dimsizes(wa) ;use dimsizes to assign the dimension sizes to individual variables
  ntim     = dsizes_x(0)
  klev     = dsizes_x(1)
  nlat     = dsizes_x(2)
  mlon     = dsizes_x(3)

;cloud level: 35:70
  lev1= 35
  lev2= 70
  lat1= 40
  lat2= 260
  lon1= 100
  lon2= 310
 
  wa_cloud= wa(:,lev1:lev2,lat1:lat2,lon1:lon2)
printMinMax(wa_cloud,1)

  ntimn = ntim
  klevn = lev2-lev1+1
  nlatn = lat2-lat1+1
  mlonn = lon2-lon1+1  
printVarSummary(wa_cloud)
  nbin  = 61
  nbin1 = nbin-1
  npts = mlonn*nlatn*klevn

;define w_bin
  w_bin= new((/nbin1,ntimn/), float)
  w_bin@_FillValue = -9999.0
  CFTD::wacftd(mlonn,nlatn,klevn,ntimn,nbin,nbin1,npts,wa_cloud,w_bin)

  w_bin = where(w_bin .eq.0,  w_bin@_FillValue, w_bin) 
printVarSummary(w_bin)
  w_bin!1 = "time"
  w_bin&time  =times 
  w_bin@description ="Fraction of vertical velocity (m/s)"

print(w_bin(:,40))

;*********************
;create plot
;*********************
  res@cnFillOn            = True            ; turn on color
  res@gsnSpreadColors     = True            ; spread out color table
  res@gsnSpreadColorStart = 2
  res@gsnSpreadColorEnd   = -3

  res@tiXAxisString  = "Time"
  res@tiYAxisString  = "vertical velocity (m/s)"

   plot  = gsn_csm_contour(wks,w_bin,res)

frame(wks)
end