load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/cd_string.ncl"
  begin

filename = "wrfout_d01_2016-10-30_00:00:00"
;
; wtype = "png"
  wtype = "pdf"
; for png --> gives "error" for pdf plot
  wtype@wkWidth  = 2500    ; produces a better looking PNG
  wtype@wkHeight = 2500
; wks = gsn_open_wks(wtype,"obdradar")

wks = gsn_open_wks(wtype,"ic_cg")
;   mpdatabaseversion= "Earth..4"

a = addfile(filename,"r")
; Set some basic resources
  res = True
  res@MainTitle                   = "ITLN"
; res@amParallelPosF   =  0.5    ; This is the right edge of the plot.
; res@amOrthogonalPosF = -0.5    ; This is the top edge of the plot.
; res@amJust           = "TopRight"
;---Top right string

;
; these don't work here
; res@gsnRightString = "gsnRightString"
; res@gsnrightStringFontHeightF = 0.015
; res@gsnRightStringFontColor = "darkorchid"
; res@RightStringParallelPosF = 1.01

; find ntimes
  times  = wrf_user_list_times(a)  ; get times in the file
  ntimes = dimsizes(times)         ; number of times in the file
  print ("ntimes = " + ntimes)
  ter = wrf_user_getvar(a,"HGT",0)
  ter = ter - 10
   printMinMax(ter,0)
  dims=dimsizes(ter)
  print("dims(0) " + dims(0))
  print("dims(1) " + dims(1))
; ascii file
a_asc = (/"output.txt"/)
;---Read the data in as a giant 1D array
    values = asciiread(a_asc,-1,"float")
    nvals  = dimsizes(values)
    print("nvals = " + nvals)
    nx=dims(1)
    ny=dims(0)
;---Calculate ntime based on the number of values we read in.
    n_l_times = nvals - ((nvals/nx)*ny)
    print("n_times = " + n_l_times)
;   ntxy  = nblk+1
    nvars  = 10
    nxy    = nx * ny
    nblk   = nxy * nvars
; create variables to hold the ascii file variables.
; 
      n_t_vars=5
      time_array     = new((/n_t_vars/),integer)
      year     = new((/n_l_times,ny,nx/),integer)
      month    = new((/n_l_times,ny,nx/),integer)
      day      = new((/n_l_times,ny,nx/),integer)
      hour     = new((/n_l_times,ny,nx/),integer)
      minute     = new((/n_l_times,ny,nx/),integer)
      lat_a    = new((/n_l_times,ny,nx/),float)
      lon_a    = new((/n_l_times,ny,nx/),float)
      ii    = new((/n_l_times,ny,nx/),integer)
      jj    = new((/n_l_times,ny,nx/),integer)
      total_light  = new((/n_l_times,ny,nx/),float)
      ic_cg  = new((/ny,nx/),float)
    do it=0,n_l_times-1
      nstart = it*(nblk+1)+1
      nend   = nstart+nblk-1
      year(it,:,:)    = toint(reshape(values(nstart:nend:nvars),(/ny,nx/)))
      month(it,:,:)    = toint(reshape(values(nstart+1:nend:nvars),(/ny,nx/)))
      day(it,:,:)    = toint(reshape(values(nstart+2:nend:nvars),(/ny,nx/)))
      hour(it,:,:)    = toint(reshape(values(nstart+3:nend:nvars),(/ny,nx/)))
      minute(it,:,:)    = toint(reshape(values(nstart+4:nend:nvars),(/ny,nx/)))
      lat_a(it,:,:)    = (reshape(values(nstart+5:nend:nvars),(/ny,nx/)))
      lon_a(it,:,:)    = (reshape(values(nstart+6:nend:nvars),(/ny,nx/)))
      ii(it,:,:)    = toint(reshape(values(nstart+7:nend:nvars),(/ny,nx/)))
      jj(it,:,:)    = toint(reshape(values(nstart+8:nend:nvars),(/ny,nx/)))
      total_light(it,:,:)    = (reshape(values(nstart+9:nend:nvars),(/ny,nx/)))
    end do
   do it = 0,n_l_times-1,1
    time=it
    print("Working on time: " + times(it) )
    res@TimeLabel = times(it)   ; Set Valid time to use on plots
;   do j=0,dims(0)-1
;   do i=0,dims(1)-1
;   print("lat_a " + lat_a(it,j,i))
;   print("lon_a " + lon_a(it,j,i))
;   end do
;   end do
    do j=0,dims(0)-1
    do i=0,dims(1)-1
    ic_cg(j,i)=total_light(it,j,i)
;   if (ic_cg(j,i).ne.0)then
;   print("light " + total_light(it,j,i))
;   end if
    end do
    end do
    max_light=max(ic_cg)
    if (max_light.eq.0)then
     ic_cg(1,1) = 1.1
    end if

  opts_ic_cg = res

  opts_ic_cg = True

;HHMM UTC D Mon YYYY
; time_array(0)=hour(1,1,1)
; time_array(1)=minute(1,1,1)
; time_array(2)=day(1,1,1)
; time_array(3)=month(1,1,1)
; time_array(4)=year(1,1,1)
; format = "HHMM UTC D Mon YYYY"
; format = ""
; time_var=cd_string(time_array,format)
; print("time_var = " + time_var)
print("year(1,1,1) = "+year(1,1,1))
print("month(1,1,1) = "+month(1,1,1))
print("day(1,1,1) = "+day(1,1,1))
print("hour(1,1,1) = "+hour(1,1,1))
print("minute(1,1,1) = "+minute(1,1,1))

new_time_units = "hours since 1800-01-01 00:00"
new_time_units = new_time_units 
new_time_array = cd_inv_calendar( year(it,1,1),  month(it,1,1), day(it,1,1), hour(it,1,1), minute(it,1,1), 0, new_time_units,0)
print("new_time_array = " + new_time_array)
printVarSummary(new_time_array)
format = "" ;; defaults to "%H%M IST %d %c %Y"
format = "%H%M IST %d %c %Y"
; 2.16667 = + 2 hours + 10 minutes (for first time on map in IST)
new_time_array = new_time_array + 2.16667
time_var=cd_string(new_time_array,format)
print("time_var = " + time_var)
  opts_ic_cg@FieldTitle           = time_var
  opts_ic_cg@gsnStringFontHeightF = 0.025
  opts_ic_cg@UnitLabel            = "10 minute ETLN Lightning Rate"
  opts_ic_cg@UnitLabelHeightF           = 0.015
  opts_ic_cg@cnLevelSelectionMode = "ExplicitLevels"
  opts_ic_cg@cnFillOn             = True


   opts_ic_cg@cnLevels             = (/ 000, 1,2., \
                                     5, 10.,15.,20,25,40,45,50/)
;  opts_ic_cg@cnFillColors         = (/"White","White","AntiqueWhite","AntiqueWhite2", \
   opts_ic_cg@cnFillColors         = (/"White","White","LightSkyBlue3","LightSkyBlue", \
                                    "SpringGreen","SpringGreen4", \
                                    "Yellow","Yellow3","Orange","Red","HotPink3","HotPink1","HotPink","Violet"/)
  opts_ic_cg@NoHeaderFooter = True
  opts_ic_cg@Footer = False



 
 contour_ic_cg = wrf_contour(a,wks,ic_cg,opts_ic_cg)
  pltres            = True   ; Set plot options
  pltres@PanelPlot  = True   ; Tells wrf_map_overlays not to remove contours
  mpres             = True   ; Set map options (none set here)
  mpres@mpOutlineOn                 = True
  mpres@mpUSStateLineColor          = "Black"
  mpres@mpNationalLineColor         = "Black"
  mpres@mpGeophysicalLineColor      = "Black"
  mpres@mpUSStateLineThicknessF     = 3.0
  mpres@mpNationalLineThicknessF    = 3.0
  mpres@mpGeophysicalLineThicknessF = 3.0
  mpres@mpOutlineBoundarySets       = "AllBoundaries"


    plot = wrf_map_overlays(a,wks,(/contour_ic_cg/),pltres,mpres)


    draw(plot)
    frame(wks)

      delete(opts_ic_cg)
   end do  ;       end of the time loop
end