; ncl 'wrfvar="PBLH"' 'plnm="pblhnew"' wrfout_ncl.ncl

;----------------------------------------------------------------------
; wrf_gsn_1.ncl
;----------------------------------------------------------------------
; Concepts illustrated:
;   - Using gsn_csm_contour_map to plot WRF-ARW data
;----------------------------------------------------------------------
; The key to making gsn_csm_contour_map work with WRF data is:
; 
; 1. Reading XLAT/XLONG off the file and attaching them as 
;    special "lat2d" and "lon2d" attributes of the data to 
;    be plotted.
;
; 2. Zooming in on the map (WRF data is regional and NCL
;    creates global plots by default).
;----------------------------------------------------------------------; 
; These files are loaded by default in NCL V6.2.0 and newer
 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
 load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
;----------------------------------------------------------------------
;main code

begin
;---Open WRF output file.
  dir      = "./"	;
  filename = systemfunc("ls -1 wrfout_d02_20*")
  print(filename)
  a = addfile(dir + filename,"r")

;---Read terrain height and lat/lon off file.
; do it=1,120,6
  it        = 120     ; first time step
;	 function wrf_user_getvar( file_handle, varin[*]:string, time:integer )
  var       = wrf_user_getvar(a,wrfvar ,it)  ; Terrain elevation
  var@lat2d = wrf_user_getvar(a,"XLAT" ,it)  ; latitude/longitude
  var@lon2d = wrf_user_getvar(a,"XLONG",it)  ; required for plotting
  printVarSummary(var)
  print ("min/max var = " + min(var) + "/" + max(var))
  
  times  = wrf_user_list_times(a)  ; get times in the file
  ntimes = dimsizes(times)         ; number of times in the file
  print("There are " +ntimes+ " times in this file")
;  print( times )
  
  wks = gsn_open_wks("png",plnm+it)

;---Set some basic plot options
  res               = True

  res@gsnMaximize   = True   ; maximize plot in frame

  res@tiMainString  = "PBLH for "+filename

  res@gsnFrame      = False
  res@cnFillOn      = True  
  res@cnFillPalette = "BlGrYeOrReVi200"
  res@lbBoxEndCapStyle = "TriangleBothEnds"
  res@cnLinesOn     = False
  res@cnLevelSelectionMode = "ExplicitLevels"
  res@cnLevels             = (/10,25,50,100,200,400,600,800,1000,1200,1400/) ; "m"
  res@lbOrientation        = "vertical"  ; vertical label bar
  res@pmLabelBarWidthF     = 0.05		; Width of labelbar
  res@lbLabelFontHeightF   = 0.012		; Size of labarbar font

;---Zoom in on plot
   minlon					 = 42.
   maxlon					 = 66.
   minlat					 = 24.
   maxlat					 = 42.
   res@mpMinLatF             = minlat ; 24. ;20. ;23. ;min(var@lat2d)
   res@mpMaxLatF             = maxlat ; 41. ;45. ;42. ;max(var@lat2d)
   res@mpMinLonF             = minlon ; 43. ;38. ;42. ;min(var@lon2d)
   res@mpMaxLonF             = maxlon ; 65. ;70. ;65. ;max(var@lon2d)
   res@mpCenterLonF          = 0.5*(res@mpMinLonF + res@mpMaxLonF)
   res@mpOutlineBoundarySets = "National"

;---Boundaries
   res@mpDataSetName            = "Earth..4"   ; 1 through 4
   res@mpNationalLineColor      = "black" ; "mpCountyLineColor" doesn't work
   res@mpNationalLineThicknessF = 2.5     ; "mpCountyLineThicknessF" doesn't work
   
;---tichmarks
   res@tmXBLabelFontHeightF     = 0.01		  ; Size of tickmarks
;   res@tmYLLabelFontHeightF     = 0.01		  ; Size of tickmarks
   res@tmXTOn             = False
   res@tmYROn             = False
   ; how set tick labels
   res@tmXBMode           = "Explicit"
;   res@tmXBMinorValues    = ispan(toint(minlon),toint(maxlon),1)
   res@tmXBMinorValues    = fspan(minlon,maxlon,23)
   lonvalues			  = fspan(minlon,maxlon,13)
   lonlabels			  = (/"42E","44E","46E","48E","50E","52E","54E","56E","58E","60E","62E","64E","66E"/)
   res@tmXBMinorOn		  = True
;   res@tmXBMinorLengthF	  = 0.5	
;   res@tmXBMajorOutwardLengthF = 0.02
;   res@tmXBMinorOutwardLengthF = 0.02
;   res@tmXBMinorLineColor = 2
;   res@tmXBMajorLineColor = 0
   res@tmXBValues         = lonvalues
   res@tmXBLabels         = lonlabels
   res@tmXBMajorLengthF   = 0.012     ; Sets the length of the X-Axis bottom major tick marks.

   res@tmYLMode           = "Explicit"
;   res@tmYLMinorValues    = fspan(minlat,maxlat,19)
   latvalues			  = fspan(minlat,maxlat,10)
   latlabels			  = (/"24N","26N","28N","30N","32N","34N","36N","38N","40N","42N"/)
   res@tmYLValues         = latvalues
   res@tmYLLabels         = latlabels
   res@tmYLMajorLengthF   = 0.012     ; Sets the length of the Y-Axis left major tick marks.

   contour = gsn_csm_contour_map(wks,var,res)
   draw(contour)
   frame(wks)
; end do

end