diri = "./"
  fili = "WRFOUT_d03_2017-06-28.T2.nc" 
  pthi = diri+fili
  a    = addfile(pthi,"r")
  
  time_dc=fspan(0,23,24)

  TIME   = a->XTIME               ; *ALL* times on the file
  nTIME  = dimsizes(TIME)         ; each time is one hour
  
  tStrt  = 72                     ; skip 0->71 [72 hrs]
  tLast  = nTIME-2                ; skip Aug 1 00Z
  print("tStrt="+tStrt+"  tLast"+tLast)

  Time   = TIME(tStrt:tLast)      ; explore time variable
  ymdh   = cd_calendar(Time,-3)   ; yyyymmddhh     
 ;print(ymdh)                            
  printMinMax(ymdh,0)             ; min=2017070100   max=2017073123

  ntimes = dimsizes(Time)
  ndays  = ntimes/24
  print("ntimes="+ntimes+"  ndays="+ndays)

  T2     = a->T2(tStrt:tLast,:,:)
  T2&Time= Time                   ; associate Time as coordinate variable
  printVarSummary(T2)
  print("------------------------")

  dimT2  = dimsizes(T2)
  ntim   = dimT2(0)               ; all hours
  nlat   = dimT2(1)
  mlon   = dimT2(2)

;---The following is the hourly climatology for the current month
  T2_dc  = T2(:23,:,:)            ; create array to hold hrly climatology
  printVarSummary(T2_dc)          ; T2_dc(24,:,:)

  do hr=0,23                      ; at each 'hr' there will be 'ndays' values 
     T2_dc(hr,:,:)= (/ dim_avg_n(T2(hr::24,:,:),0) /)
  end do
  
  T2_dc@long_name = "Hourly climatology for "+(ymdh(0)/10000)
  printVarSummary(T2_dc)          ; [Time | 24] x [south_north | 303] x [west_east | 297]
  printMinMax(T2_dc,0)
  print("------------------------")
  print("num(ismissing(T2_dc)="+num(ismissing(T2_dc))) ; gives the number of missing values in T2_dc
  print("------------------------")

;---Computing the composite diurnal cycle of anomalies
  T2mean = dim_avg_n_Wrap(T2_dc,0)  ; ???????????????  ; effectively the is is the monthly mean
  T2mean@long_name = "T2mean = dim_avg_n_Wrap(T2_dc,0)  <== MONTHLY MEAN"
  printVarSummary(T2mean)         ; [south_north | 303] x [west_east | 297]
  printMinMax(T2mean,0)           ; min=284.294   max=299.132
  print("------------------------")

;======???=====> Do not understand why the following two lines
;;T2_dc_anom = T2_dc
;;T2_dc_anom = (/ T2_dc-conform(T2_dc,T2mean,(/1,2/)) /)  ; difference from overall monthly mean

  T2_dc_anom = T2      ; create array with meta data
  T2_dc_anom@long_name = "Anomalies from Hourly Climatology"

  do hr=0,23           
     T2_work = T2(hr::24,:,:)     ; [Time | 31] x [south_north | 303] x [west_east | 297]
     T2_dc_anom(hr::24,:,:) = T2_work - conform(T2_work, T2_dc(hr,:,:), (/1,2/))
  end do 
  T2_dc_anom@long_name = "Hourly anomalies: "+(ymdh(0)/10000)

  printVarSummary(T2_dc_anom)     ; [Time | 744] x [south_north | 303] x [west_east | 297]
  printMinMax(T2_dc_anom,0)       ; min=-7.73807   max=6.89316      
  print("------------------------")
  delete(T2_work)