;   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"

begin
;
; We generate plots, but what kind do we prefer?
  type = "x11"
; type = "pdf"
; type = "ps"
; type = "ncgm"

; We know the data ia big, so we will increase the workstation space
  wks = gsn_open_wks(type,"olr")
  setvalues NhlGetWorkspaceObjectId()
  "wsMaximumSize": 33554432
  end setvalues

; Set some Basic Plot Information
; We would like to switch off the Initial Time and Footer info on the plot
; We also do not want titles on the color label bar
  ARWres = True
  ARWres@MainTitle                   = "NRCM WRF"
  ARWres@InitTime                    = False
  ARWres@Footer                      = False
  ARWres@lbTitleOn                   = False

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

  a = addfile("wrfout_d02_2005-08-26_00:00:00.nc","r")

  times  = wrf_user_list_times(a)  ; get times in the file
  ntimes = dimsizes(times)         ; number of times in the file
  print("Working on OLR")

; Keep this for Valid time on the plot
  ARWres@TimeLabel = times(0)

; We are interested in a zoomed area with the SW corner at 20N;85W,  
; and the NE corner at 37N;40W
; So get the XY points of these points
  lats = (/  20.0,  37.0 /)
  lons = (/ -85.0, -40.0 /)
  loc = wrf_user_latlon_to_ij(a, lats, lons)
 
; loc(;,0) is south-north (y) and loc(:,1) is west-east (x)
  y_start = loc(0,0)
  y_end = loc(1,0)
  x_start = loc(0,1)
  x_end = loc(1,1)


; Change color map to something that has a grey scale
  gsn_define_colormap(wks,"wxpEnIR")

; Create a basic and zoom map background
  mpres = True
  map = wrf_map(wks,a,mpres)                   
  map_zoom = wrf_map_zoom(wks,a,mpres,y_start,y_end,x_start,x_end)

  olr = wrf_user_getvar(a,"OLR",0)
  olr_zoom = olr(y_start:y_end,x_start:x_end)

  opts = ARWres                         
  opts@FieldTitle = "OLR"
  opts@cnFillOn = True  
  opts@ContourParameters = (/ 100., 340., 20./)
  opts@gsnSpreadColorStart = 51 
  opts@cnRasterModeOn = True
  contour  = wrf_contour(a,wks,olr,opts)
  contourZ = wrf_contour(a,wks,olr_zoom,opts)

; Plot OLR in grey scale for the full domain
  wrf_map_overlay(wks,map,(/contour/),True)

; Do it again, but do not advance the frame, so we can plot
; a box on the output
  wrf_map_overlay(wks,map,(/contour/),False)
        lnres = True
        lnres@gsLineColor = "NavyBlue"
        lnres@gsLineThicknessF = 3.0

        gsn_polyline(wks,map,(/-85.0,-40.0/),(/20.0,20.0/),lnres)
        gsn_polyline(wks,map,(/-85.0,-40.0/),(/37.0,37.0/),lnres)
        gsn_polyline(wks,map,(/-85.0,-85.0/),(/20.0,37.0/),lnres)
        gsn_polyline(wks,map,(/-40.0,-40.0/),(/20.0,37.0/),lnres)

        frame(wks)  ; Now that we are done drawing, draw the frame

; Now zoom into the box area and plot OLR just for this area
  wrf_map_overlay(wks,map_zoom,(/contourZ/),True)
  delete(opts)
  delete(mpres)
  delete(map_zoom)

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

  print("Working on RAIN ")

; Open a second input file 
  b = addfile("auxhist1_d02_2005-08-25_01:00:00.nc","r")

; Change color map
  gsn_define_colormap(wks,"WhViBlGrYeOrReWh")

; Change some of the default map background options
  mpres = True
  mpres@mpGeophysicalLineColor      = "Brown"
  mpres@mpNationalLineColor         = "Brown"
  mpres@mpUSStateLineColor          = "Brown"
  mpres@mpGeophysicalLineThicknessF = 2.0
  mpres@mpNationalLineThicknessF    = 2.0
  mpres@mpUSStateLineThicknessF     = 2.0
  map_zoom = wrf_map_zoom(wks,a,mpres,y_start,y_end,x_start,x_end)

; Get surface wind and slp
  u10 = wrf_user_getvar(a,"U10",0)
  v10 = wrf_user_getvar(a,"V10",0)
    u10_zoom = u10(y_start:y_end,x_start:x_end)
    v10_zoom = v10(y_start:y_end,x_start:x_end)
  slp = wrf_user_getvar(a,"slp",0)
    wrf_smooth_2d( slp, 3 )
    slp_zoom = slp(y_start:y_end,x_start:x_end)

; Get rainfall for two times
  rainc_25 = wrf_user_getvar(b,"RAINC",22)
  rainnc_25 = wrf_user_getvar(b,"RAINNC",22)
  rainc_26 = wrf_user_getvar(b,"RAINC",23)
  rainnc_26 = wrf_user_getvar(b,"RAINNC",23)

; Calculate rainfall tendencies between these two times
  rrn = (rainc_26+rainnc_26) - (rainc_25+rainnc_25)
  rrn_zoom = rrn(y_start:y_end,x_start:x_end)


  opts = ARWres                         
  opts@FieldTitle = "Hourly Acc Rainfall"
  opts@cnLevelSelectionMode = "ExplicitLevels"
  opts@cnLevels = (/ .1, .2, .4, .8, 1.6, 3.2, 6.4, \
                      12.8, 25.6, 51.2, 102.4/)
  opts@cnFillColors = (/"White","White","DarkOliveGreen1", \
                        "DarkOliveGreen3","Chartreuse", \
                        "Chartreuse3","Green","ForestGreen", \
                        "Yellow","Orange","Red","Violet"/)
  opts@cnInfoLabelOn = False
  opts@cnConstFLabelOn = False
  opts@cnFillOn = True  
  opts@cnRasterModeOn = True
  opts@gsnSpreadColorEnd  = -3 
  cont_rain = wrf_contour(a,wks,rrn_zoom,opts)
  delete(opts)

  opts = ARWres
  opts@ContourParameters = (/ 980., 1024., 4. /)
  opts@cnLineColor = "NavyBlue"
  opts@cnLineLabelsOn = False
  opts@gsnContourLineThicknessesScale = 2.0
  cont_slp = wrf_contour(a,wks,slp_zoom,opts)
  delete(opts)

  opts = ARWres
  opts@FieldTitle = "Wind"
  opts@NumVectors = 47
  vector = wrf_vector(a,wks,u10_zoom,v10_zoom,opts)
  delete(opts)

; Plot rainfall and wind vectors and then
; rainfall, wind vectors and slp
  wrf_map_overlay(wks,map_zoom,(/cont_rain,vector/),True)
  wrf_map_overlay(wks,map_zoom,(/cont_rain,cont_slp,vector/),True)

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

end