diri    = "./"                             ; input directory   
   fili    = "MCD19A2.A2008*.hdf"                         ; input file       
   f       = addfile (diri+fili , "r")                          
;***********************************************************
; Read time and create required yyyyddd                    
;***********************************************************
   time    = f->time                          ; time:units = "hours since 1-1-1 00:00:0.0"                               
   TIME    = cd_calendar(time, 0)          ; type float 
   year    = toint( TIME(:,0) )
   month   = toint( TIME(:,1) )
   day     = toint( TIME(:,2) ) 
   ddd     = day_of_year(year, month, day) 
   yyyyddd = year*1000 + ddd                                  ; needed for input
;***********************************************************
; Read data: short2flt                                     
;*********************************************************** 
   aod     = short2flt( f->Optical_Depth_047_grid1km(:,0,:,:) )    ; convert to float 
   printVarSummary( aod ) 
;***********************************************************
; Compute daily climatology: raw and then 'smoothed'  
;***********************************************************
   hClmDay    = clmDayTLL(aod, yyyyddd)     ; daily climatology at each grid point                                      
   printVarSummary(hClmDay)                                                                                              
;***********************************************************
; Compute smoothed daily climatology using 2 harmonics
;***********************************************************
   hClmDay_sm = smthClmDayTLL(hClmDay, 2)
   printVarSummary(hClmDay_sm)

;***********************************************************
; Compute daily anomalies using raw and smoothed daily climatologies
;***********************************************************
    hAnom      = calcDayAnomTLL (aod, yyyyddd, hClmDay)     
    printVarSummary(hAnom_sm)
    printMinMax(hAnom, 0)

    hAnom_sm   = calcDayAnomTLL (aod, yyyyddd, hClmDay_sm)     

    hAnom_sm@long_name = "Anomalies from Smooth Daily Climatology"
    printVarSummary(hAnom_sm)
    printMinMax(hAnom_sm, 0)