;----------------------------------------------------------------------
; wrf_ozone.ncl
;----------------------------------------------------------------------
; Concepts illustrated:
;   - Using wrf_xxxx scripts to plot WRF-ARW data
;----------------------------------------------------------------------
; This script is meant to show the difference between plotting WRF
; data using wrf_xxx scripts, and using gsn_csm_xxx scripts.
;
; See wrf_gsn_4.ncl for an example of using gsn_csm_contour_map to
; plot WRF data.
;----------------------------------------------------------------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

begin
;---Open WRF output file
;   filename = "wrfout_d01_2010-06-30_12:00:00.nc"
;   filename = "wrfout_d01_2010-06-30_13:00:00.nc"
;   filename = "wrfout_d01_2010-06-30_14:00:00.nc"
;   filename = "wrfout_d01_2010-06-30_15:00:00.nc"
;   filename = "wrfout_d01_2010-06-30_16:00:00.nc"
;   filename = "wrfout_d01_2010-06-30_17:00:00.nc"
;   filename = "wrfout_d01_2010-06-30_18:00:00.nc"
;   filename = "wrfout_d01_2010-06-30_19:00:00.nc"
;   filename = "wrfout_d01_2010-06-30_20:00:00.nc"
;   filename = "wrfout_d01_2010-06-30_21:00:00.nc"
;   filename = "wrfout_d01_2010-06-30_22:00:00.nc"
;   filename = "wrfout_d01_2010-06-30_23:00:00.nc"
;   filename = "wrfout_d01_2010-07-01_07:00:00.nc"
;   filename = "wrfout_d01_2010-07-01_08:00:00.nc"
    filename = "wrfout_d01_2010-07-05_20:00:00.nc"
;  filename = "dailyave_2007-06-01.nc"
   
  a = addfile(filename,"r")

  
; Set some Basic Plot options
  res                       = True
  res@MainTitle             = "REAL-TIME WRF"
  res@lbLabelAutoStride     = True            ; let NCL determine label spacing
  res@mpDataBaseVersion     = "MediumRes"
  res@mpDataSetName         = "Earth..4"      ; U.S. counties
  InitTitle                 = False
  
;  res@lbOrientation        = "Horizontal"
;  lbTitlePosition = "Right"

  mpres = True
  mpres@mpGeophysicalLineColor    = "Black"
  mpres@mpNationalLineColor       = "Blue"
  mpres@mpUSStateLineColor        = "Black"
  mpres@mpCountyLineColor         = "grey50"
  mpres@mpDataBaseVersion         = "MediumRes"
  mpres@mpOutlineBoundarySets     = "AllBoundaries"
  mpres@mpUSStateLineThicknessF   = 2.0
  mpres@mpNationalLineThicknessF  = 2.5
  mpres@mpDataResolution          ="Finest"
  mpres@mpProjection              ="LambertConformal"
  mpres@pmTickMarkDisplayMode     = "Always"
  
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; What times and how many time steps are in the data set?
  times = wrf_user_getvar(a,"times",-1)  ; get all times in the file
  ntimes = dimsizes(times)        ; number of times in the file

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

  do it = 00,ntimes -1,1           ; TIME LOOP
    
    print("Working on time: " + times(it) )
    res@TimeLabel = times(it)   ; Set Valid time to use on plots

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

; First get the variables we will need        
;---Read temperature at first time step

    o3  = wrf_user_getvar(a,"o3",it)    ; O3
    u   = wrf_user_getvar(a,"ua",it)    ; 3D U at mass points
    v   = wrf_user_getvar(a,"va",it)    ; 3D V at mass points
    
    nl  = 0 ; bottomost level
;    o3 = o3(nl,:,:)
    u10 = u(nl,:,:)
    v10 = v(nl,:,:)
    
    o3=o3*1000.
    o3@description = "Surface Ozone Concentration"
    o3@units="ppbv"
    u10 = u10*1                            ; Keeps wind units in m/s
    v10 = v10*1
    u10@units       = "m/s"
    v10@units       = "m/s"

   wks = gsn_open_wks("pdf","ozone_emiss")

;---Resources for filled contour plot
  
;   res@ContourParameters = (/0., 68., 4./)    ; Change the spacing
;   res@cnLevelSelectionMode="AutomaticLevels"
;   res@cnMinLevelValF = 0
;   res@cnMaxLevelValF = 50
;   res@cnLevelSpacingF = 2

   nl   = 0  ; bottomost level
   ntimes = 12 ; number of times to be accounted for in script
  
  ;---Wind vector plot
  vec_res            = res         
;  vec_res@FieldTitle = "x-wind component"       ; overwrite Field Title
  vec_res@NumVectors = 30         ; density of wind barbs
  vec_res@vcLineArrowThicknessF = 0.1
  vec_res@vcGlyphStyle          = "CurlyVector"

  vector = wrf_vector(a,wks,u10,v10,vec_res)
  
;  plot = wrf_contour(a,wks,o3(nl,:,:),res)


; 03
  opts = res
  opts@lbLabelAutoStride   = True            ; let NCL determine label spacing
  opts@lbLabelBarOn        = True
  opts@lbTitlePosition     = "Bottom"
  opts@cnLineColor         = "Blue"
  opts@ContourParameters   = (/0.,68.,4./) ; Change the spacing
  opts@cnFillOn            = True
  opts@gsnContourLineThicknessesScale = 0.5
  contour = wrf_contour(a,wks,o3(nl,:,:),opts)
  plot = wrf_map_overlays(a,wks,(/contour,vector/),res,mpres)
  delete(opts)

;---Overlay plot on map and draw.
;  ov = wrf_map_overlays(a,wks,plot,True,True)

 end do

end