external GET_IJ "./get_ij.so"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
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/contrib/cd_string.ncl"
 begin
 tel_aviv = True
 plots_collated = new(4,graphic)

;station_num,date,time (LST),station_name,longitude,latitude,height (m),tas,rh,wind_spd (m/s),wind_dir,p_station
delim = ","
  wrf_file ="../RUN/wrfout_d03_2018-07-23_00:00:00_NEW_LAI_1.5_SMCDRY"
  a = addfile(wrf_file,"r")
  xlat2d = wrf_user_getvar(a,"XLAT",0)
  xlon2d = wrf_user_getvar(a,"XLONG",0)
  xland = wrf_user_getvar(a,"XLAND",0)
  lu_index = wrf_user_getvar(a,"LU_INDEX",0)
  times  = wrf_user_list_times(a)  ; get times in the file
  ntimes = dimsizes(times)         ; number of times in the file
  print ("ntimes = " + ntimes)

  t2 = wrf_user_getvar(a,"T2",-1)
  t2@_FillValue = -999
  printVarSummary(t2)
; print(x)
  rh2 = wrf_user_getvar(a,"rh2",-1)
; wspd_wdir10; 6.6 or higher
  u10 = wrf_user_getvar(a,"U10",-1)
  u10 = wrf_user_getvar(a,"U10",-1)
  v10 = wrf_user_getvar(a,"V10",-1)
  speed = sqrt(u10^2+v10^2)
  wdir = wind_direction(u10,v10,0)
  heat_index = wrf_user_getvar(a,"AFWA_HEATIDX",-1)
  dims2d=dimsizes(xlat2d)
  asci_file = "IMS_10_Min_data_JER.csv"
; asci_file = "IMS_10_Min_data_TA.csv"
  lines = asciiread(asci_file,-1,"string")
  csv_lines := str_split_csv(lines(1:),delim,0)
  nlines = dimsizes(lines)-1
  print("nlines = " + nlines)
  Name = new(nlines,string)
  Latitude = new(nlines,float)
  Longitude = new(nlines,float)
  Date = new(nlines,string)
  Hour = new(nlines,integer)
  Minute = new(nlines,integer)
  Day = new(nlines,integer)
  Year= new(nlines,integer)
  Month= new(nlines,integer)
  Height = new(nlines,integer)
  Temp = new(nlines,float)
  Heat = new(nlines,float)
  Wind_Speed = new(nlines,float)
  Wind_Dir = new(nlines,float)
  Pressure = new(nlines,float)
  Rh = new(nlines,float)
  time_var = new(ntimes,string)
  t_o = new(ntimes,float)
  t_m = new(ntimes,float)
  hi_o = new(ntimes,float)
  hi_m = new(ntimes,float)
  ws_o = new(ntimes,float)
  ws_m = new(ntimes,float)
  wd_o = new(ntimes,float)
  wd_m = new(ntimes,float)
; 23,23/07/2018,0:00,JERUSALEM CENTRE,35.2217,31.7806,810,22.3,60,3.9,278,917.5
; do i_line = 0,nlines-1
    Station                         := csv_lines(:,0)    ; read column 1
  Date               = csv_lines(:,1)    ; read column 2
; print("Date" + Date)
  Year= tointeger(str_get_cols(Date, 6, 9))
; print(Year)
  Month= tointeger(str_get_cols(Date, 3,4))
  Day= tointeger(str_get_cols(Date, 0, 1))
  TimeLST            = csv_lines(:,2)    ; read column 2
; print("TimeLST  = " + TimeLST)
  Hour= tointeger(str_get_cols(TimeLST, 00, 01))
  Minute= tointeger(str_get_cols(TimeLST, 03, 04))
; print(Year)
; print(Month)
; print(Day)
; print(Hour)
  Longitude            = tofloat(csv_lines(:,4))    ; read column 2
; print(Longitude)
  Latitude             = tofloat(csv_lines(:,5))    ; read column 2
; printVarSummary(Latitude)
; print(Latitude)
  Height             = tointeger(csv_lines(:,6))    ; read column 2
; print(Height)
  Temp             = tofloat(csv_lines(:,7))    ; read column 2
; print(Temp)
; print(x)
  Rh             = tofloat(csv_lines(:,8))    ; read column 2
; print(Rh)
  Wind_Speed             = tofloat(csv_lines(:,9))    ; read column 2
; print(Wind_Speed)
  Wind_Dir             = tofloat(csv_lines(:,10))    ; read column 2
; print(Wind_Dir)
  Pressure             = tofloat(csv_lines(:,11))    ; read column 2
  io = (/2,2/)                                       ; degF input, degF output
  t = Temp*1.8+32
  rh = Rh
  H_I = heat_index_nws(t, rh, io, False)
  H_I = (H_I-32)/1.8
; printMinMax(H_I,False)


; print(Pressure)
  ;end do
   ii_loc = 0
   jj_loc = 0
   GET_IJ::get_ij(xlat2d,xlon2d,Latitude(0),Longitude(0),ii_loc,jj_loc,dims2d(0),dims2d(1))
;  ii_loc = ii_loc-1
   ii_loc = ii_loc+3
   jj_loc = jj_loc-1
   print("ii_loc = " + ii_loc)
   print("jj_loc = " + jj_loc)
   print(xland(jj_loc,ii_loc))
;  print(x)
   it_val = 0
   Second = Year
   Second = 0
   i_val = 0
   printMinMax(heat_index,False)
   heat_index = where(heat_index.eq.0,avg(heat_index),heat_index)
   printMinMax(heat_index,False)
;where(lu_index.eq.31.or.lu_index.eq.32.or.lu_index.eq.33,tgr_1_ave,tgr_null)
   t2 = t2 - 273.13
   heat_index = heat_index - 273.13
   dims2d = dimsizes(lu_index)
      
   do it_val = 0,nlines,6
   it = i_val
   tgr_null = t2(it,:,:)
   tgr_null = t2@_FillValue
   if (tel_aviv.eq.True)then
   t2_lu=where((lu_index.eq.31.or.lu_index.eq.32.or.lu_index.eq.33).and.xlon2d.gt.35,t2(it,:,:),tgr_null)
   hi_lu=where((lu_index.eq.31.or.lu_index.eq.32.or.lu_index.eq.33).and.xlon2d.gt.35,heat_index(it,:,:),tgr_null)
   ws_lu=where((lu_index.eq.31.or.lu_index.eq.32.or.lu_index.eq.33).and.xlon2d.gt.35,speed(it,:,:),tgr_null)
   wd_lu=where((lu_index.eq.31.or.lu_index.eq.32.or.lu_index.eq.33).and.xlon2d.gt.35,wdir(it,:,:),tgr_null)
   end if
    t_m(i_val) = avg(t2_lu)
    t_o(i_val) = Temp(it_val)
;   hi_m(i_val) = heat_index(i_val,jj_loc,ii_loc) - 273.13
    hi_m(i_val) = avg(hi_lu)
;   printMinMax(heat_index(i_val,jj_loc,ii_loc),False)
;   print(x)
    hi_o(i_val) = H_I(it_val)
;   ws_m(i_val) = speed(i_val,jj_loc,ii_loc)
    ws_m(i_val) = avg(ws_lu)
    ws_o(i_val) = Wind_Speed(it_val)
;   wd_m(i_val) = wdir(i_val,jj_loc,ii_loc)
    wd_m(i_val) = avg(wd_lu)
    wd_o(i_val) = Wind_Dir(it_val)
    new_time_units = "hours since 1800-01-01 00:00"
    new_time_units = new_time_units
    new_time_array = cd_inv_calendar( Year(it_val),  Month(it_val), Day(it_val), Hour(it_val), Minute(it_val), Second(it_val), new_time_units,0)
    print("minute = " + Minute(it_val))
    print("second = " + Second(it_val))
;   print("new_time_array = " + new_time_array)
    format = "" ;; defaults to "%H%M EST %d %c %Y"
;   format = "" ;; defaults to "%H%M EST %d %c %Y"
;   format = "%H%M%S IST %D %N %Y"
    format = "%H IST %D %N"
    new_time_array = new_time_array + 0
    print("new_time_array = " + new_time_array)
    time_var(i_val)=cd_string(new_time_array,format)
    print("time_var = " + time_var(i_val))
    print("t_m = " + t_m(i_val))
    print("t_o = " + t_o(i_val))
    print("HI_m = " + hi_m(i_val))
    print("HI_o = " + hi_o(i_val))
    print("ws_m = " + ws_m(i_val))
    print("ws_o = " + ws_o(i_val))
    print("wd_m = " + wd_m(i_val))
    print("wd_o = " + wd_o(i_val))
    print("i_val = " + i_val)
    i_val = i_val + 1
   end do
   printVarSummary(heat_index)
   printMinMax(H_I,False)
;  print(x)
;---To plot multiple lines, you must put them into a mulidimensional array.
;data      = new((/3,dimsizes(variable_1(:,0))/),float)
 data      = new((/2,dimsizes(t_m)/),float)
;clim = new(4,float)
;clim(0) = 15.47
;clim(1) = 10.91
;clim(2) = 6.12
;clim(3) = 2.80
 wks   = gsn_open_wks ("png","OBS_vs_Mod_Jer")
;wks   = gsn_open_wks ("png","OBS_vs_Mod_TA")
 do i_plot = 0,3,1
 print("i_plot = " + i_plot)
 if (i_plot.eq.0)then
 data(0,:) = t_m
 data(1,:) = t_o
 end if
 if (i_plot.eq.1)then
 data(0,:) = hi_m
 data(1,:) = hi_o
 end if
 if (i_plot.eq.2)then
 data(0,:) = ws_m
 data(1,:) = ws_o
 end if
 if (i_plot.eq.3)then
 data(0,:) = wd_m
 data(1,:) = wd_o
 end if
;print("data(1,:) = " + data(0,:))
;print("data(2,:) = " + data(1,:))
 xaxis = new(dimsizes(data(0,:)-1),float)
;xaxis(0) = 0
 xaxis(:) = ispan(0,ntimes-1,1)


;---Set plotting parameters

 res                   = True                      ; plot mods desired

;res@tmXTOn      = False
 res@tmYLFormat  = "f"      ; remove trailing ".0"

;--------------------------------------------------
; The time_axis_label function adds additional
; resources to "res" to produce nicely-formatted
; time labels on X axis. This function only works
; if you have a time "units" recognized by the
; cd_calendar function.
;--------------------------------------------------
  restick           = True
  res = True
; restick@ttmFormat = "%N/%D %H:%M"
; restick@ttmFormat = "%N/%D"

; time_axis_labels(times_c,res,restick)
;;x_axis_labels(times_c,res,restick)


  n_t = dimsizes(data(0,:)) 
; res@tiMainString      = "Disruptive Snowstorms"       ; add title
  res@tiMainString      = " "
; set min/max
  min_h = min(data(:,:))
  max_h = max(data(:,:))
  res@trYMinF     = min_h-1
  res@trYMaxF     = max_h + 1
; res@trYMinF     = tointeger(min_h) - 50 		; set minimum Y-axis value
; res@trYMaxF     = tointeger(max_h) + 50  ; maximum Y-axis value
  res@trXMinF     = 20

; res@trXMaxF     =  ntimes
  res@tmXBMode = "Explicit"
  res@tmXBValues            =  xaxis(0::24)              ; Indicate where we want tickmar
  print("res@tmXBValues = " + res@tmXBValues)
; dimsx = dimsizes(x)
; xx = x(0:dimsx(0)-1)
; res@tmXBValues = time_var(0::6)
  res@tmXBLabels = time_var(0::24)
  print("res@tmXBValues = " + res@tmXBValues)
; res@tmXBLabelFontHeightF = 0.019
; print("xaxis = " + xaxis)
  res@tmXBLabelAngleF = 67.5
; print(x)
  label = (/"m","o"," Inter","Refer."," Satur."," Leafy"/)
  colors      = (/"violet","red","orange","green","darkgreen","blue"/)
  res@xyLabelMode            = "Custom"           ; label a line
  res@xyExplicitLabels       = label            ; text to use
  res@xyLineLabelFontHeightF = 0.0200              ; font height
  res@xyLineColors   = colors              ; label color
  res@xyLineLabelFontColors   = colors              ; label color
  if (i_plot.eq.0)then
  res@tiMainString         = "Tel-Aviv"
  res@tiYAxisString = "Temp (2 m)" +   " [~S~o~N~C]"
  end if
  if (i_plot.eq.1)then
  res@tiMainString         = "Tel-Aviv"
  res@tiYAxisString = "Heat Index (2 m)" +   " [~S~o~N~C]"
  end if
  if (i_plot.eq.2)then
  res@tiMainString         = "Tel-Aviv"
  res@tiYAxisString = "Wind Speed" +   " [m s~S~-1~N~]"
  end if
  if (i_plot.eq.3)then
  res@tiMainString         = "Tel-Aviv"
  res@tiYAxisString = "Wind Direction" +   " [~S~o~N~]"
  end if


;
; Similiar resources are xyLineThicknessF and xyLineColor,
; which will affect all lines in the array.
;
 res@xyLineThicknesses = (/  3.0,   3.0,3.0,3.0,3.0,3.0/)
 res@xyLineColors      = colors

 plot  = gsn_csm_xy (wks,xaxis,data,res) ; create plot
 plots_collated(i_plot) = plot
 end do
  resP                = True                ; panel only resources
  resP@gsnMaximize    = True                ; maximize plots
  resP@gsnPanelFigureStrings = (/"(e)","(f)","(g)","(h)"/)
  resP@gsnPanelFigureStringsFontHeightF = 0.02


 gsn_panel(wks,plots_collated,(/2,2/),resP)          ; now draw as one plot
    
  end