;*************************************************
; WRF: Panel an individual variable at 1 hour time steps 
;************************************************
load "/ptmp/lamptey/wrf_v2.1.2/process/WRF_NCL/WRFOptions.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/contributed.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "/ptmp/lamptey/wrf_v2.1.2/process/WRF_NCL/WRFPlot.ncl"
load "/ptmp/lamptey/wrf_v2.1.2/process/WRF_NCL/WRFUserARW.ncl"
load "/ptmp/lamptey/wrf_v2.1.2/process/WRF_NCL/SkewTFunc.ncl"
begin
;************************************************
; get all desired files
;************************************************
  diri = "/ptmp/lamptey/wrf_v2.1.2/process/data/explicit/"
  fili = systemfunc("cd "+diri +"; ls explB_thomE_out_d02*")
  fili = fili +".nc"                           ; explicitly add file extension

  F    = addfiles (diri+fili, "r")             ; refers to *all* files      
  f    = addfile  (diri+fili(1), "r")          ; refers to 2nd file only
  print(fili)
  print(F)
;************************************************
; open file and read in data
; (1) Read RAINC 
; (2) Read character variable Times; Convert to string for plots
;************************************************
  newl    = inttochar(10)    ; new line character. use to span
                             ; multiple lines for nice presentation
; x       = addfiles_GetVar (F, diri+fili,"RAINC"); read RAINC@ across all files
; x       = addfiles_GetVar (F, diri+fili,"RAINNC"); read RAINC@ across all files
; x@description = "WRF model precip at 4 km resolution" + newl +\ 
;                   "(NCEP & Kain_Fritch conv. & Lin et al. microphys.)"
;                   "(NCEP & Kain_Fritch conv. & thompson et al. microphys.)"
;                   "(NCEP & Kain_Fritch conv. & wsm6 microphys.)"
; times   = chartostring(F[:]->Times)          ; built-in function

  x = f->RAINC + f->RAINNC 
  x@description = "Hourly precip at 4 km resolution" 
  times   = chartostring(f->Times)          ; built-in function
  ntim    = dimsizes(times)                    ; # of time steps
  print(ntim)
; print(times)
  min_x = min(x)
  max_x = max(x)
; print(min_x)
; print(max_x)
;************************************************
; create plots  
;************************************************
  wks = gsn_open_wks("X11" ,"mcs_d2")           ; ps,pdf,x11,ncgm,eps
  gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ; select color map

  res                       = True             ; plot mods desired
  res@gsnSpreadColors       = True             ; use full range of colormap
; res@gsnSpreadColorStart   = 47             ; use full range of colormap
; res@gsnSpreadColorEnd     = 159             ; use full range of colormap
  res@cnFillOn              = True             ; color plot desired
  res@cnLinesOn             = False            ; turn off contour lines
  res@cnLineLabelsOn        = False            ; turn off contour labels

;************************************************
; Use WRF_contributed procedure to set map resources
;************************************************
  WRF_map_c(f, res, 0)                   ; reads info from file
  res@mpProjection      = "Mercator"       ; choose projection
  res@mpGridAndLimbOn   = False              ; turn on lat/lon lines
  res@mpPerimOn         = False             ; turn off box around plot
  res@mpGridLatSpacingF = 5.               ; spacing for lat lines
  res@mpGridLonSpacingF = 5.               ; spacing for lon lines
  res@mpFillOn          = False
  res@mpOutlineBoundarySets       = "National"   ; turn on country boundaries
  res@mpGeophysicalLineColor      = "Navy"       ; color of cont. outlines
  res@mpGeophysicalLineThicknessF = 1.5          ; thickness of outlines
  
;************************************************
; set True for native projection (faster)
;************************************************
  res@tfDoNDCOverlay       = True

;************************************************
; For individual plots 
;************************************************
  res@gsnDraw               = False            ; do not draw 
  res@gsnFrame              = False            ; do not advance 'frame'
  res@lbLabelBarOn          = False            ; turn off individual lb's 
;************************************************
; common label bar for panel requires that
; all plots should be set to the same interval
; Use built in function to determine "nice" limits
;************************************************
  mnmxint                   = nice_mnmxintvl( min(x), max(x), 14, False)
  res@cnLevelSelectionMode  = "ManualLevels"   
  res@cnMinLevelValF        = mnmxint(0)
  res@cnMaxLevelValF        = mnmxint(1)
  res@cnLevelSpacingF       = mnmxint(2) ;2
;************************************************
; Allocate array to store plots: specify time step
;************************************************
  tstep                     = 1                ; time step to plot          
  plts = new (ntim/tstep , "graphic")          ; 1d array to hold plots

;************************************************
; loop over each forecast 
;************************************************
  n    = -1                                    ; counter
  do nt=0,ntim-1,tstep                         ; plot every 1 hour
     res@gsnLeftString         = times(nt)
     n  = n+1
     plts(n)                   = gsn_csm_contour_map(wks,x(nt,:,:),res)
  end do
;************************************************
; create panel: panel plots have their own set of resources
;************************************************
  resP                  = True                   ; modify the panel plot
  resP@gsnFrame         = False
  resP@txString         = x@description 
  resP@gsnPanelLabelBar = True                   ; add common colorbar
  resP@gsnMaximize      = True                   ; maximize panel area
  resP@lbLabelAutoStride= True                   ; let NCL figure lb stride

  resP@gsnPanelTop      = 0.95                   ; create some space at top
  resP@gsnPanelBottom   = 0.10                   ; add space at bottom
; resP@gsnPanelFigureStrings= (/"a)","b)","c)",\
;                               "d)","e)","f)"/) ; add strings to panel
; resP@amJust           = "BottomLeft"           ; move labels
  gsn_panel(wks,plts,(/6,4/),resP)               ; now draw as one plot
  
  add_T               = "NCEP, BMJ (coarse domain only),., and Thompson @4 km resolution"
  xpos                = 0.5
  ypos                = 0.05

  txres               = True ; text mods desired
  txres@txFontHeightF = 0.01  ; choose text font
  gsn_text_ndc(wks,add_T,xpos,ypos,txres)
  frame(wks)
end