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/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "/data/muhdomer/NSCC/scratch/WRF_myfunc.ncl"

begin

wrf_file1 = systemfunc("ls /data/muhdomer/Build_WRF/Martilli/glade/u/home/martilli/wrf_3.8/WRFNCEP_updatedUFRAC/wrfout_d05_2016-04-01_00:00:00")
  wrf_file2 = systemfunc("ls  /data/muhdomer/Build_WRF/Martilli/glade/u/home/martilli/wrf_3.8/WRF_Forest/run/wrfout_d05_2016-04-01_00:00:00")

  a = addfiles(wrf_file1+".nc","r")
  b = addfiles(wrf_file2+".nc","r")
  ListSetType(a,"cat")
  ListSetType(b,"cat")

  LH1 = addfiles_GetVar(a,wrf_file1,"T2")
  LH2 = addfiles_GetVar(b,wrf_file2,"T2")
  
  LH = LH1 - LH2
  ;LH=LH-273.16 
  
  times = a[:]->Times
  ntimes = dimsizes(times)         ; number of times in the file
  ntimes_plot = ntimes(0)  


      
  print(ntimes_plot)
  lh3d = new((/29,129,210/), float)
  lh2d = new((/129,210/), float)
    c = 0 
    do i=time_serie, ntimes_plot-2,24
       lh2d(:,:) = LH(i,:,:)
       lh3d(c,:,:) = lh2d(:,:)
       c = c + 1
           
      end do
    lh_avg = dim_avg_n(lh3d,0)

  

  lh_avg@description = "UHI"
  
  
  lh_avg@units = "C"
  ; Plotting options for LH
      type = "png"
      wks = gsn_open_wks(type,"LH_d05_diff_3_1_"+hour)
      res =True  
      opts = res
      opts@cnFillOn = True
     
      opts@cnLevelSelectionMode = "ExplicitLevels"
       opts@cnLevels=(/-0.2,0,0.2,0.4,0.6,0.8,1,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,4.0,4.2,4.4,4.6/)
     
      pltres = True
      pltres@NoTitles = True
      pltres@gsnMaximize = True
      pltres@FramePlot = False
      pltres@PanelPlot = True
      mpres = True
      contour_lh = wrf_contour(a[0],wks,lh_avg,opts)
      delete(opts)
   
   
    u10 = addfiles_GetVar(a,wrf_file1,"U10")
    v10 = addfiles_GetVar(a,wrf_file1,"V10")
    
    u103d = new((/29,129,210/), float)
    u102d = new((/129,210/), float)
    c = 0
    do i=time_serie+6, ntimes_plot-1,24
       
       u102d(:,:) = u10(i,:,:)
       u103d(c,:,:) = u102d(:,:)
       c = c + 1

     end do
    u10_avg = dim_avg_n(u103d,0)
   

    v103d = new((/29,129,210/), float)
    v102d = new((/129,210/), float)
    c = 0
    do i=time_serie+6, ntimes_plot-1,24
       v102d(:,:) = v10(i,:,:)
       v103d(c,:,:) = v102d(:,:)
       c = c + 1

    end do
    v10_avg = dim_avg_n(v103d,0)
 
    
     opts=res
     opts@vcGlyphStyle               = "LineArrow"
     opts@vcLineArrowThicknessF      = 3
     vector  = wrf_vector(a[0],wks,u10_avg,v10_avg,opts)
     delete(opts)
       
      

 plot = wrf_map_overlays(a[0],wks,(/contour_lh,vector/),pltres,mpres)


   shape_files = (/"/data/muhdomer/NSCC/scratch/Shapefiles/SGP_adm0.shp","/data/muhdomer/NSCC/scratch/Shapefiles/MYS_adm0.shp","/data/muhdomer/NSCC/scratch/Shapefiles/IDN_adm0.shp"/)

    plot = add_coastline_sgmyid(shape_files,wks,plot,"black")

     draw(plot)
     frame(wks)
end