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"

begin
  a = addfile("/data/muhdomer/Build_WRF/Martilli/glade/u/home/martilli/wrf_3.8/WRFNCEP_updatedUFRAC/wrfout_d05_2016-04-01_00:00:00.nc","r")   ; the file used to extract index information for each land use type
  lu = wrf_user_getvar(a, "LU_INDEX", 0)   ; the land use index
  lu_oned = ndtooned(lu)                   ; one dimension version of land use index
  dims = dimsizes(lu)

  index_CHighR = ind_resolve(ind(lu_oned .eq. 31),dims)   
  index_CMidR = ind_resolve(ind(lu_oned .eq. 32), dims) 
  index_CLowR = ind_resolve(ind(lu_oned .eq. 33), dims) 
  index_OHighR = ind_resolve(ind(lu_oned .eq. 34),dims)        
  index_OMidR = ind_resolve(ind(lu_oned .eq. 35), dims) 
  index_OLowR = ind_resolve(ind(lu_oned .eq. 36), dims) 
  index_LWLowR = ind_resolve(ind(lu_oned .eq. 37),dims)
  index_LLowR = ind_resolve(ind(lu_oned .eq. 38), dims) 
  index_SB = ind_resolve(ind(lu_oned .eq. 39), dims) 
  index_HI = ind_resolve(ind(lu_oned .eq. 40), dims)


  index_forest = ind_resolve(ind(lu_oned .eq. 2), dims) ; evergreen broadleaf forest
  index_shrub = ind_resolve(ind(lu_oned .eq. 7), dims)  ; open shrubland
  index_grass = ind_resolve(ind(lu_oned .eq. 10), dims) ; grasslands
  index_barren = ind_resolve(ind(lu_oned .eq. 16), dims) ; barren or sparsely vegetated

  dims_CHighR = dimsizes(index_CHighR)
  dims_CMidR = dimsizes(index_CMidR)
  dims_CLowR = dimsizes(index_CLowR)
  dims_OHighR = dimsizes(index_OHighR)
  dims_OMidR = dimsizes(index_OMidR)
  dims_OLowR = dimsizes(index_OLowR)
  dims_LWLowR = dimsizes(index_LWLowR)
  dims_LLowR = dimsizes(index_LLowR)
  dims_SB = dimsizes(index_SB)
  dims_HI = dimsizes(index_HI)



  dims_forest = dimsizes(index_forest)
  dims_shrub = dimsizes(index_shrub)
  dims_grass = dimsizes(index_grass)
  dims_barren = dimsizes(index_barren)


  type = "x11"
; type = "pdf"
; type = "eps"
;type = "newpng"
  wks = gsn_open_wks(type,"UHI_2016_April")



  times  = wrf_user_list_times(a)  ; get times in the file
  ntimes = dimsizes(times)         ; number of times in the file

  temp = new((/14,24/),float)      ; average surface temperature for each land use type with diurnal cycle
;
;  do it = 0,23             ; TIME LOOP
  do it = 1, ntimes-1,23
;
   print("Working on time: " + it )


    tc2 = wrf_user_getvar(a,"T2",it)     ; T2 in Kelvin
    tc2 = tc2-273.16                     ; T2 in C

    it_local = (it+8)%24
; print("Local time: " + it_local )

    tc2_CHighR= 0.0
    do i=0, dims_CHighR(0)-1
       tc2_CHighR = tc2_CHighR + tc2(index_CHighR(i,0), index_CHighR(i,1))

    end do

    temp(0, it_local) = tc2_CHighR/dims_CHighR(0)


    tc2_CMidR = 0.0
    do i=0, dims_CMidR(0)-1
       tc2_CMidR = tc2_CMidR + tc2(index_CMidR(i,0), index_CMidR(i,1))
    end do

    temp(1, it_local) = tc2_CMidR/dims_CMidR(0)


    tc2_CLowR = 0.0
    do i=0, dims_CLowR(0)-1
       tc2_CLowR = tc2_CLowR + tc2(index_CLowR(i,0), index_CLowR(i,1))
    end do
   print("it_local="+it_local+" forest= "+temp(2,it_local))
    temp(2, it_local) = tc2_CLowR/dims_CLowR(0)

 tc2_OHighR = 0.0
    do i=0, dims_OHighR(0)-1
       tc2_OHighR = tc2_OHighR + tc2(index_OHighR(i,0), index_OHighR(i,1))

    end do

    temp(3, it_local) = tc2_OHighR/dims_OHighR(0)

tc2_OMidR = 0.0
    do i=0, dims_OMidR(0)-1
       tc2_OMidR = tc2_OMidR + tc2(index_OMidR(i,0), index_OMidR(i,1))

    end do

    temp(4, it_local) = tc2_OMidR/dims_OMidR(0)

tc2_OLowR = 0.0
    do i=0, dims_OLowR(0)-1
       tc2_OLowR = tc2_OLowR + tc2(index_OLowR(i,0), index_OLowR(i,1))

    end do

    temp(5, it_local) = tc2_OLowR/dims_OLowR(0)

tc2_LWLowR = 0.0
    do i=0, dims_LWLowR(0)-1
       tc2_LWLowR = tc2_LWLowR + tc2(index_LWLowR(i,0), index_LWLowR(i,1))

    end do

    temp(6, it_local) = tc2_LWLowR/dims_LWLowR(0)

tc2_LLowR = 0.0
    do i=0, dims_LLowR(0)-1
       tc2_LLowR = tc2_LLowR + tc2(index_LLowR(i,0), index_LLowR(i,1))

    end do

    temp(7, it_local) = tc2_LLowR/dims_LLowR(0)

tc2_SB = 0.0
    do i=0, dims_SB(0)-1
       tc2_SB = tc2_SB + tc2(index_SB(i,0), index_SB(i,1))

    end do

    temp(8, it_local) = tc2_SB/dims_SB(0)

tc2_HI = 0.0
    do i=0, dims_HI(0)-1
       tc2_HI = tc2_HI + tc2(index_HI(i,0), index_HI(i,1))

    end do

    temp(9, it_local) = tc2_HI/dims_HI(0)




    tc2_lu_forest = 0.0
    do i=0, dims_forest(0)-1
       tc2_lu_forest = tc2_lu_forest + tc2(index_forest(i,0), index_forest(i,1))
    end do
;    tc2_lu_forest = tc2_lu_forest/dims_forest(0)
    temp(10, it_local) = tc2_lu_forest/dims_forest(0)
    print("it_local="+it_local+" forest= "+temp(10,it_local))

    tc2_lu_shrub = 0.0
    do i=0, dims_shrub(0)-1
       tc2_lu_shrub = tc2_lu_shrub + tc2(index_shrub(i,0), index_shrub(i,1))
    end do
;    tc2_lu_shrub = tc2_lu_shrub/dims_shrub(0)
    temp(11, it_local) = tc2_lu_shrub/dims_shrub(0)
    print("it_local="+it_local+" shrub= "+temp(11,it_local))

    tc2_lu_grass = 0.0
    do i=0, dims_grass(0)-1
       tc2_lu_grass = tc2_lu_grass + tc2(index_grass(i,0), index_grass(i,1))
    end do
;    tc2_lu_grass = tc2_lu_grass/dims_grass(0)
    temp(12, it_local) = tc2_lu_grass/dims_grass(0)
   print("it_local="+it_local+" grass= "+temp(12,it_local))

    tc2_lu_barren = 0.0
    do i=0, dims_barren(0)-1
       tc2_lu_barren = tc2_lu_barren + tc2(index_barren(i,0), index_barren(i,1))
    end do
;    tc2_lu_barren = tc2_lu_barren/dims_barren(0)
    temp(13, it_local) = tc2_lu_barren/dims_barren(0)
    print("it_local="+it_local+" barren= "+temp(13,it_local))

  end do        ; END OF TIME LOOP

    uhi_intensity = new((/10, 24/), float)
    uhi_intensity(0, :) = temp(0, :) - temp(11,:)
    uhi_intensity(1, :) = temp(1, :) - temp(11,:)
    uhi_intensity(2, :) = temp(2, :) - temp(11,:)
    uhi_intensity(3, :) = temp(3, :) - temp(11,:)
    uhi_intensity(4, :) = temp(4, :) - temp(11,:)
    uhi_intensity(5, :) = temp(5, :) - temp(11,:)
    uhi_intensity(6, :) = temp(6, :) - temp(11,:)
    uhi_intensity(7, :) = temp(7, :) - temp(11,:)
    uhi_intensity(8, :) = temp(8, :) - temp(11,:)
    uhi_intensity(9, :) = temp(9, :) - temp(11,:)





    ; MAKE PLOTS                                       
    x = ispan(0, 23, 1)
    res = True
    res@gsnMaximize = True
 
    res@gsnDraw              = False
    res@gsnFrame             = False

    res@tmXBLabelsOn = True
    
    res@xyLineThicknesses = (/2.0, 2.0, 3.0/)
;    res@xyLineColors    =  (/"red", "blue", "black"/)
    res@xyLineColors    =  (/"red", "blue", "black", "green", "cyan", "purple", "saddlebrown","violetred4", "violet", "turquoise"/)
    res@tiXAxisString = "Local time (h)"
    res@tiYAxisString = "UHI intensity (~S~o~N~C)"
;
;    ; add a legend
    res@pmLegendDisplayMode    = "Always"              ; turn on legend

    res@pmLegendSide           = "Top"                 ; Change location of
    res@pmLegendParallelPosF   = 0.55                 ; move units right
    res@pmLegendOrthogonalPosF = -0.4                  ; more neg = down

    res@pmLegendWidthF         = 0.2                  ; Change width and
    res@pmLegendHeightF        = 0.2                  ; height of legend.
    res@lgLabelFontHeightF     = .010                ; change font height
    res@lgPerimOn              = False                 ; no box around
    
 
    res@xyExplicitLegendLabels = (/"CompactHR","CompactMR","CompatLR","OpenHR","OpenMR","OpenLR","LWLR","LLR","SB","HI"/)
    res@lgItemOrder  = (/9,8,7,6,5,4,3,2,1,0 /)   ; Reorder the legends
;
    plot = gsn_csm_xy(wks, x, uhi_intensity, res)
;    plot = gsn_csm_xy(wks, x, temp, res)

    gres = True
    plot_line = gsn_add_polyline(wks, plot, (/0.0, 24.0/), (/0.0, 0.0/), gres)
    draw(plot)
    frame(wks)
end