;   Example script to produce plots for a WRF real-data run,
;   with the ARW coordinate dynamics option.
;   Script show how to zoom into a given area

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

begin
;
; The WRF ARW input file.  
; This needs to have a ".nc" appended, so just do it.
  a = addfile("../wrfout_d01_2000-01-24_12:00:00.nc","r")

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

; Set some Basic Plot options
  ARWres = True
  ARWres@MainTitle                   = "REAL-TIME WRF"

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

; What times and how many time steps are in the data set?
  FirstTime = True
  times  = wrf_user_list_times(a)  ; get times in the file
  ntimes = dimsizes(times)         ; number of times in the file

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

  do it = 0,ntimes-1,2             ; TIME LOOP

    print("Working on time: " + times(it) )
    ARWres@TimeLabel = times(it)   ; Set Valid time to use on plots

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  
    tc  = wrf_user_getvar(a,"tc",it)     ; T in C
    ter = wrf_user_getvar(a,"HGT",it)    ; terrain height

    dims = dimsizes(ter)
    ; As an example, we are looking for the lower right 1/4 of the domain
    x_start = dims(1)/2
    x_end   = dims(1)-1
    y_start = 0
    y_end   = dims(0)/2

    tc_zoom  = tc(:,y_start:y_end,x_start:x_end)
    tc_zoom@description = tc@description + " - Lower right of domain"
    ter_zoom = ter(y_start:y_end,x_start:x_end)
    ter_zoom@description = ter@description + " - Lower right of domain"

    mpres = True                   ; Create map backgrounds
    map      = wrf_map(wks,a,mpres) 
    map_zoom = wrf_map_zoom(wks,a,mpres,y_start,y_end,x_start,x_end)

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

    if ( FirstTime )
      ; Plot full terrain, then zoomed terrain
      opts = ARWres                    
      opts@cnFillOn = True  
      contour = wrf_contour(a,wks,ter,opts)
      wrf_map_overlay(wks,map,(/contour/),True)
      contour = wrf_contour(a,wks,ter_zoom,opts)
      wrf_map_overlay(wks,map_zoom,(/contour/),True)
      delete(opts)
      FirstTime = False
    end if


    ; Plot full tc, then zoomed tc
      opts = ARWres                    
      opts@cnFillOn = True  
      contour = wrf_contour(a,wks,tc(0,:,:),opts)
      wrf_map_overlay(wks,map,(/contour/),True)
      contour = wrf_contour(a,wks,tc_zoom(0,:,:),opts)
      wrf_map_overlay(wks,map_zoom,(/contour/),True)
      delete(opts)

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

  end do        ; END OF TIME LOOP

end