;   Example script to produce an OLR plot for a WRF real-data run,
;   with the ARW coordinate dynamics option.
;   In this example we are also going to zoom into an area of interest

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

begin
;   
; The WRF ARW input files  
; This needs to have a ".nc" appended, so just do it.
  a = addfile("./wrfout_d01_2005-08-28_00.nc","r")


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

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


; Set some Basic Plot Information
  res = True
  res@MainTitle    = "KATRINA"   ; name of plot
  res@InitTime     = False       ; do not plot initial time
  res@Footer       = False	 ; switch footers off

; Set up two sets of plot and map resources
  pltres = True
  pltres1 = True
  pltres1@FramePlot = False      ; For this resorce - do not advance frame automatically
  mpres = True
  mpres1 = True

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

  times = wrf_user_getvar(a,"times",-1)  ; get all times in the file
  ntimes = dimsizes(times)         ; number of times in the file

; We are interested in a zoomed area with the SW corner at 20N;90W,  
; and the NE corner at 30N;80W
; So get the XY points of these points
  lats = (/  20.0,  30.0 /)
  lons = (/ -93.0, -80.0 /)
  loc = wrf_user_ll_to_ij(a, lons, lats, True)
 
; loc(0,;) is west-east (x) ; loc(1,:) is south-north (y) 
; subtract one since we want to use it as an index in NCL
  x_start = loc(0,0) - 1
  x_end   = loc(0,1) - 1
  y_start = loc(1,0) - 1
  y_end   = loc(1,1) - 1

  mpres1@ZoomIn = True        ; set up map info for zoomed area
  mpres1@Xstart = x_start
  mpres1@Ystart = y_start
  mpres1@Xend = x_end
  mpres1@Yend = y_end


  olr = wrf_user_getvar(a,"OLR",-1)              ; get olr
  olr_zoom = olr(:,y_start:y_end,x_start:x_end)  ; create a zoomed area

  opts = res                         
  opts@FieldTitle = "OLR"                        ; overwrite the filed name
  opts@cnFillOn = True  
  opts@lbTitleOn = False                         ; remove field name from label bar
  opts@ContourParameters = (/ 100., 340., 20./)
  opts@gsnSpreadColorStart = 51 

  do it = 1,ntimes-1                             ; start at time 2

    opts@TimeLabel = times(it)                   ; set Valid time on the plot
  
    contour  = wrf_contour(a,wks,olr(it,:,:),opts)
    contourZ = wrf_contour(a,wks,olr_zoom(it,:,:),opts)

; Plot OLR in grey scale for the full domain
    plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)

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

        gsn_polyline(wks,plot,(/lons(0),lons(1)/),(/lats(0),lats(0)/),lnres)
        gsn_polyline(wks,plot,(/lons(0),lons(1)/),(/lats(1),lats(1)/),lnres)
        gsn_polyline(wks,plot,(/lons(0),lons(0)/),(/lats(0),lats(1)/),lnres)
        gsn_polyline(wks,plot,(/lons(1),lons(1)/),(/lats(0),lats(1)/),lnres)

        frame(wks)  ; Now that we are done with the box, advance the frame

; Now zoom into the box area and plot OLR just for this area
    plot = wrf_map_overlays(a,wks,(/contourZ/),pltres,mpres1)

  end do

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

end