;   Example script to produce plots for a WRF real-data run,
;   with the ARW coordinate dynamics option.
;   The script shows 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"

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
  res = True
  res@MainTitle                   = "REAL-TIME WRF"

  pltres = True
  mpres = True
  mpres1 = True


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

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

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

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

    print("Working on time: " + times(it) )
    res@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

    mpres1@ZoomIn = True
    mpres1@Xstart = x_start
    mpres1@Ystart = y_start
    mpres1@Xend = x_end
    mpres1@Yend = y_end

    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"

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

    if ( FirstTime )
      ; Plot full terrain, then zoomed terrain
      opts = res                    
      opts@cnFillOn = True  
      contour = wrf_contour(a,wks,ter,opts)
      plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)
      contour = wrf_contour(a,wks,ter_zoom,opts)
      plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres1)
      delete(opts)
      FirstTime = False
    end if


    ; Plot full tc, then zoomed tc
      opts = res                    
      opts@cnFillOn = True  
      contour = wrf_contour(a,wks,tc(0,:,:),opts)
      plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)
      contour = wrf_contour(a,wks,tc_zoom(0,:,:),opts)
      plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres1)
      delete(opts)

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

  end do        ; END OF TIME LOOP

end