begin

;================================================
; read data
;================================================

  f = addfile("data/netcdf/clo20210209.nc","r")

  time   = f->epochTime
  height = f->height
  bcksct = f->Backscatter

  printVarSummary(time)
  printVarSummary(height)
  printVarSummary(bcksct)

  ndims = dimsizes(bcksct)
  nt = ndims(0)
  print("nt = "+nt)

  inv_bcksct = bcksct(height|:,time|:)
  printVarSummary(inv_bcksct)

;============================================================
;-- convert the time proleptic_gregorian calendar to UTC date

  year    = new((/nt/),"integer")
  month   = new((/nt/),"integer")
  day     = new((/nt/),"integer")
  hour    = new((/nt/),"integer")
  minute  = new((/nt/),"integer")
  second  = new((/nt/),"integer")

  julianday  = new((/nt/),"integer")
  doy        = new((/nt/),"integer")

  date_str = new((/dimsizes(time)/),"string")
  date_str_i = new((/dimsizes(time)/),"integer")

  do i=0,nt-1

    utc_date = cd_calendar(time(i), 0)

;-- set date variable names

    year(i)      = tointeger(utc_date(:,0))
    month(i)     = tointeger(utc_date(:,1))
    day(i)       = tointeger(utc_date(:,2))
    hour(i)      = tointeger(utc_date(:,3))
    minute(i)    = tointeger(utc_date(:,4))
    second(i)    = tointeger(utc_date(:,5))

;   date_str(i) = sprinti("%0.4i",year(i))   + "-" + \
;                 sprinti("%0.2i",month(i))  + "-" + \
;                 sprinti("%0.2i",day(i))    + "-" + \
;                 sprinti("%0.2i",hour(i))   + "-" +  \
;                 sprinti("%0.2i",minute(i)) + "-" + \
;                 sprinti("%0.2i",second(i))

    date_str(i) = sprinti("%0.2i",hour(i))   + ":" + \
                  sprinti("%0.2i",minute(i)) + ":" + \
                  sprinti("%0.2i",second(i))

    if(i.le.30) then
      print("date_str("+i+")= "+date_str(i))
    end if
 
  end do

  printVarSummary(date_str)

;============================================================
; plot
;============================================================

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

  res                 = True                ; plot mods desired

  res@gsnMaximize     = True

  res@vpWidthF        = 0.8                 ; set width and height
  res@vpHeightF       = 0.4

  res@sfXArray        = date_str
; res@sfXArray        = time
; res@sfYArray        = height

  res@cnInfoLabelOn   = False               ; turn off contour info label

  res@cnFillMode      = "RasterFill"
  res@cnRasterSmoothingOn  = True

  res@cnFillOn        = True                ; turn on color
  res@cnFillPalette   = "gui_default"       ; set color map
  res@cnFillPalette   = "BlAqGrYeOrReVi200"       ; set color map
  res@cnLineLabelsOn  = False               ; turn off line labels
  res@cnLinesOn       = False               ; turn off contour lines

  res@tiMainString    = "Ceilometer backscatter"
  res@tiMainFontHeightF  = 0.02 

; --- X axis ---

  res@tiXAxisString            = "Time"
  res@tiXAxisFontHeightF       = 0.014 

  res@tmXBLabelFontHeightF     = 0.012
  res@tmXBMajorLengthF         = 0.012
  res@tmXBLabelAngleF          = 90.

; --- Y axis ---

  res@tiYAxisString            = "Height (m)"
  res@tiYAxisFontHeightF       = 0.014 

  res@tmYLLabelFontHeightF     = 0.012 
  res@tmYLMajorLengthF         = 0.012

  res@pmLabelBarWidthF          = 0.7
  res@pmLabelBarHeightF         = 0.1
  res@lbLabelFontHeightF        = 0.012

  plot  = gsn_csm_contour(wks,inv_bcksct,res)         ; contour the variable

end