load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"      ; These two libraries are automatically
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"    ; loaded from NCL V6.4.0 onward.
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl"


;----------------------------------------------------------------------
; This function combines two date arrays that are assumed to have
; the same units and calendar.
;----------------------------------------------------------------------
function combine_dates(date1,date2)

begin
;
; This code makes sure you span the full range of both
; dates and that you cover any gaps in the middle.
;
  min_date = min((/min(date1),min(date2)/))
  max_date = max((/max(date1),max(date2)/))
  new_date = ispan(toint(min_date),toint(max_date),1)
  new_date@units    = date1@units
  new_date@calendar = date1@calendar
  return(new_date)
end


begin

;1)******************** READ CSV  and SET TIME ******************************************

;1.a)----------------> PLUME HEIGHT 

 plume_file = asciiread("Plume_height_ECV.dat", -1, "string")  ; -1 means read whole file

;DEBUG
;print(plume_file)                               ; echo input
;exit 
 
 str     = str_split_csv (plume_file, ",", 0)
 str@_FillValue = -999.0                  ; change to 'nicer' _FillValue

 plume_H = tofloat(str(:,1))
 plume_H@_FillValue = -999                  ; change to 'nicer' _FillValue

;DEBUG
;print(str)
;print(plume_H)
;exit 
 
 delim = ","
 nfields = str_fields_count(plume_file(0), delim)   ; Count the fields separated
                                                    ; by one or more spaces.
;DEBUG
;print(nfields)                                     ; nfields = 2
;exit

 field = 1
 subStrings = str_get_field(plume_file, field, delim)

;DEBUG
;print(subStrings)
;exit

 hour = toint(str_get_cols(subStrings, 0, 1))

;DEBUG
;print(hour)
;exit

 minute = toint(str_get_cols(subStrings, 3, 4))

;DEBUG
;print(minute)
;exit

 second = toint(str_get_cols(subStrings, 6, 7))

;DEBUG
;print(second)
;exit

 ndates = dimsizes(hour)
 year   = new(ndates,integer)
 year   = 2013
 month = new(ndates,integer)
 month = 11
 day   = new(ndates,integer)
 day   = 23


;DEBUG
;print("=============================")
;print(year+"_"+month+"_"+day+"_"+hour+"_"+minute+"_"+second)
;exit

;1.b)----------------> MER file 

 MER_file = asciiread("MER.dat", -1, "string")  ; -1 means read whole file

;DEBUG
;print(MER_file)                               ; echo input
;exit 
 
 str2    = str_split_csv (MER_file, ",", 0)
 str2@_FillValue = 0                  ; change to 'nicer' _FillValue

 MER_value = tofloat(str2(:,1))
 MER_value@_FillValue = 0                  ; change to 'nicer' _FillValue
 ;MER_value = where(MER_value.eq.0, MER_value@_FillValue, MER_value)

;DEBUG
;print(str2)
;print(MER_value)
;exit 
 
 delim2 = ","
 nfields2 = str_fields_count(MER_file(0), delim2)   ; Count the fields separated
                                                    ; by one or more spaces.
;DEBUG
;print(nfields2)                                     ; nfields = 2
;exit

 field2 = 1
 subStrings2 = str_get_field(MER_file, field2, delim2)

;DEBUG
;print(subStrings2)
;exit

 hour2= toint(str_get_cols(subStrings2, 0, 1))

;DEBUG
;print(hour2)
;exit

 minute2= toint(str_get_cols(subStrings2, 3, 4))

;DEBUG
;print(minute2)
;exit

 second2= toint(str_get_cols(subStrings2, 6, 7))

;DEBUG
;print(second2)
;exit

 ndates2= dimsizes(hour2)
 year2   = new(ndates2,integer)
 year2   = 2013
 month2 = new(ndates2,integer)
 month2 = 11
 day2   = new(ndates2,integer)
 day2   = 23

;print("=============================")
;print(year2+"_"+month2+"_"+day2+"_"+hour2+"_"+minute2+"_"+second2)
;exit

;2) **********************************************************************
;Convert dates on CSV file to same units/calendar as WRF file
 opt = 1
 units  = "hours since 1900-01-01 00:00:00" ; "seconds/hours/days since ...."
 csv_time  = cd_inv_calendar (year,month,day,hour,minute,second,units,opt)
 csv_time2 = cd_inv_calendar(year2,month2,day2,hour2,minute2,second2,units,opt)

 ;---Combine the two date arrays for a single date array for the X axis labels.
 axis_date = combine_dates(csv_time,csv_time2)


;DEBUG
;print(axis_date)
;print(csv_time)
;print(csv_time2)
;exit

;3) ************************************************************************
; Workstation settings

  type ="png"
  wks = gsn_open_wks(type,"fig4") ; send graphics to PNG file

  res             = True
;  res@gsnMaximize = True
  res@gsnDraw     = False
  res@gsnFrame    = False

  res@vpWidthF    = 0.8      ; Make plots wider than 
  res@vpHeightF   = 0.2      ; they are high.

  res@tmYLOn               = True     ; Turn off left tickmarks
  res@tmYROn               = False    ; Turn on right tickmarks
  res@tmXTOn               = False   ; Turn off top tickmarks

  res@tmYLLabelsOn         = True    ; Turn off left labels
  res@tmYRLabelsOn         = True    ; Turn on right labels
  res@tmYRMinorOn          = False   ; Turn off minor ticks on Y axis
  res@tmXBMinorOn          = False   ; Turn off minor ticks on Y axis

  res@tmYLMinorOn          = True 

  res@tmYLLabelFontHeightF = 0.01   ; Increase font height
;  res@tmYLLabelDeltaF      = 2.0     ; Increase space b/w ticks and labels
  res@tmYRLabelJust        = "CenterRight"   ; right-justify labels
  
;  res@trXMinF          = min(csv_time)
;  res@trXMaxF          = 998387

;-------- Surface

  res@tiXAxisString = "Time [h]"
  res@xyMarkers        = 4           ; marker styles
  restick           = True
  restick@ttmFormat = "%H:%M:%S"
  restick@ttmMajorStride =  30  
  restick@ttmMinorStride = 1
  time_axis_labels(csv_time,res,restick)
  print(restick)

; Change y axis string and color for each plot.
  res@tiYAxisString = "xy1"
  res@xyLineThicknessF = 5
  res@xyLineColor   = "blue"


; Overrule time_axis_labels settings
  print(res)   
  ln = 998383   
  csv_time0 = (csv_time-ln)
  res@trXMinF          = min(csv_time) - ln
  res@trXMaxF          = 998387-ln
  res@tmXBMinorValues = res@tmXBMinorValues-ln
  res@tmXBValues = res@tmXBValues - ln
  res@tmXBOn = False
  res@tmXBLabelsOn = False
  res@tmXBMinorValues := fspan(1.25,4.75,8)
  res@tmXBMinorOn = True
  res@trYMinF = -0.5
  print(res)

  plot0 = gsn_csm_xy(wks,csv_time0,plume_H,res)

  res@tmXBOn = True
  res@tmXBLabelsOn = True
;  res@trYLog = True
 res@xyYStyle          = "Log"
 res@tmYLMode          = "Explicit"                ; explicit labels
 res@tmYLValues        = (/10,1000,10000,100000,10^6,10^7/)
 res@tmYLLabels        = ""+res@tmYLValues 
 res@xyMarkLineMode   = "Markers"          ; markers and lines
 res@tiYAxisString = "xy2"

 res@trYMinF = 10
 res@trYMaxF = max(res@tmYLValues)
 csv_time3 = (csv_time2-ln)
 print(res)
  plot1  = gsn_csm_xy(wks,csv_time3,MER_value,res)

; Set up resource lists for attaching the plot.
; The res1 will apply to the base plot, and the
; res2 to the plots being attached.  These resources
; lists are *not* for changing things like line color,
; but for changing things like whether the plots
; are maximized, and which axis they are attached on.
;
  res1                     = True
  res2                     = True
  res1@gsnMaximize         = True
  res2@gsnAttachPlotsXAxis = True

; xy1 will be the base plot.
  amid = gsn_attach_plots(plot0,plot1,res1,res2)

  draw(plot0)    ; All four plots will be drawn.
  frame(wks)


end