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/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
begin
  ;; read data
  ;; systemfunc: Executes a shell command and returns the output.
  ;; Retrieve a list of files from one directory
  wrfout = systemfunc("ls /public3/jtwiles/Tubbs_Fire_Analysis/TUBBS_FIRE_WRFOUT/wrfout_d01_2017-10-09_05:00:00")  
  filnum = dimsizes(wrfout) 
  data   = addfiles(wrfout, "r") 
  wrflat = data[0]->XLAT    
  wrflon = data[0]->XLONG
;  printVarSummary(wrflon)
;  printVarSummary(wrflat) 
  dims   = dimsizes(wrflat)
  ;; get the dimension size of wrflat
  nt     = dims(0)
  nlat   = dims(1) 
  nlon   = dims(2)
 
  res=True
  res@gsnDraw=False               ;-- don't draw plot yet
  res@gsnFrame=False              ;-- don't advance frame
  res@gsnMaximize=False           ;-- don't maximized the workstation
  res@gsnLeftString=""
  res@gsnRightString=""
  res@mpOutlineBoundarySets="GeophysicalAndUSStates" 
  res@mpOutlineOn                 = True   ;-- outline land (default: False)
  res@mpGeophysicalLineThicknessF = 2      ;-- thickness of continental
  res@mpNationalLineThicknessF    = 2      ;-- thickness of national boundary outlines
  res@mpFillOn                    = False  ;-- no map fill 
  res@mpProjection           = "LambertConformal"  ;-- Transform using the conical Lambert Conformal projection.
  res@mpLambertParallel1F    = 43.427
  res@mpLambertParallel2F    = 43.427
  res@mpLambertMeridianF     = -112.185
  res@mpUSStateLineColor = "gray35"       ;-- color of interior provincial or state boundary outlines
  res@mpNationalLineColor = "gray35"      ;-- color of interior national boundary outlines.
  res@mpGeophysicalLineColor = "gray35"   ;-- color of all continents, islands, and inland water bodies
  res@mpLimitMode="Corners"
  res@mpLeftCornerLatF=wrflat(0,5,5)              ;-- min latitude
  res@mpLeftCornerLonF=wrflon(0,5,5)              ;-- min longitude
  res@mpRightCornerLatF=wrflat(0,nlat-5,nlon-5)   ;-- max latitude
  res@mpRightCornerLonF=wrflon(0,nlat-5,nlon-5)   ;-- max longitude

  res@pmTickMarkDisplayMode="Always"              ;--always displays TickMark object
  res@tmXBLabelFontHeightF=0.016                  ;-- sets font height of bottom X Axis labels
  res@tmYLLabelFontHeightF=0.016                  ;-- sets font height of Y-Axis left labels
  res@tmXTOn=False    ;-- turns off top tick marks
  res@tmYROn=False    ;-- turns off right tick marks
  res@gsnAddCyclic = False
  res@tfDoNDCOverlay = True
  res@cnInfoLabelOn  = False               ;-- not draw an informational label
  res@cnLinesOn            = True          ;-- draw contour lines
  res@cnLineLabelsOn       = True          ;-- draw contour lines labels
  res@cnLineLabelBackgroundColor = "Transparent"       ;-- draw lines background color
  res@cnLineColor = "black"                ;-- color for all contour lines
  res@cnLineThicknessF = 2                 ;-- thickness for all contour lines
  res@cnLineDashPattern = 0                ;-- dash pattern for all contour lines.0: solid lines

  res@cnLevelSelectionMode = "ManualLevels"  ;-- set contour levels manually
  res@cnMinLevelValF = -10                   ;-- minimum contour level 
  res@cnMaxLevelValF = 10                  ;-- maximum contour level
  res@cnLevelSpacingF = 1                    ;-- contour level spacing
  res@cnLineLabelFontHeightF = 0.016         ;-- set font height of line labels. 
  res@cnLineLabelFontColor = "black"         ;-- set color for all contour line labels.
  res@cnLinesOn      = False
  res@cnFillOn = True
  res@cnFillPalette = "BlWhRe"
;  res@vcLevelPalette = ""

  res@cnLineLabelPlacementMode = "constant"   ;-- draws line labels as an integral part of the line dash pattern
  res@cnLineLabelInterval  = 2   ;-- set label intervals on contours
  ; res@cnLineLabelDensityF  = 1   ;-- set number of labels along a contour line
  res@cnLineDashSegLenF    = 0.2

  do i = 0, filnum-1
      ; Get the data
    p     = wrf_user_getvar(data[i],"pressure", -1)
    tk = wrf_user_getvar(data[i],"tk", -1)
    u = wrf_user_getvar(data[i],"ua", -1)
    v = wrf_user_getvar(data[i],"va",-1)
    Time   = wrf_user_getvar(data[i], "times", -1) ;-- get variable : Times in file (return strings)
    print(Time)
   ; printVarSummary(T)

;   Interpolate variables to desired level
    T   = wrf_user_intrp3d(tk,p,"h",850.,0.,False)
    U   = wrf_user_intrp3d(u,p,"h",850,0.,False)
    V   = wrf_user_intrp3d(v,p,"h",850,0.,False)

  
printVarSummary(T)


  ntim  = 1
;  klev  = 1
  lat   = wrflat(0,:,0)
  lon   = wrflon(0,0,:)
;   Calculate dT/dX in Temp. Advection Equation

   dlon = (lon(2)-lon(1))*0.0174533 ; convert to radians
                                    ; pre-allocate space
   dTdX = new ( (/ntim,nlat,nlon/), typeof(T), T@_FillValue)
 printVarSummary(dTdX)
;;   do ni=0,nlon-1                      ; loop over each longitude
     do nj=0,nlat-1				; loop over each latitude

      dX = 6378388.*cos(0.0174533*lat(nj))*dlon  ; constant at this latitude
      dTdX(:,nj,:) = center_finite_diff(T(:,nj,:) , dX , False,0)
   end do
;;  end do
;  result: dTdX(time,lev,lat,lon)


;  Calculate dT/dY in Temp. Advection Equation

  dlat = (lat(2)-lat(1))*0.0174533 ; convert to radians
                                    ; pre-allocate space
   dTdY = new ( (/ntim,nlat,nlon/), typeof(T), T@_FillValue)

   do ni=0,nlon-1                      ; loop over each longitude
;;    do nj=0,nlat-1
      dY = 6378388.*cos(0.0174533*lat)*dlat  ; array, vary at different latitude
      dTdY(:,:,ni) = center_finite_diff(T(:,:,ni) , dY , False,0)
;;   end do
  end do
; result: dTdY(time,lev,lat,lon)

;  Calculate Temperature Advection

   Tadvect = (-1*(U*dTdX + V*dTdY))

     do j = 0, dimsizes(Time)-1
      wks=gsn_open_wks("x11", "Tadvect"+Time) 
    res@gsnLeftString  = "Temp. Advect. 850(hPa)"        
    res@gsnRightString  = Time(j)
      plot=gsn_csm_contour_map(wks, Tadvect(j,:,:), res) ;-- plot data using gsn_csm_contour_map function 
           draw(plot)
      frame(wks)
     end do
      delete(Tadvect)
      delete(p)
      delete(Time)
      delete(tk)
      delete(u)
      delete(v)
      delete(T)
      delete(U)
      delete(V)
   end do

end