;*******************************************
; crcm_1.ncl
;
; Concepts illustrated:
;   - Reading multiple CRCM netCDF files containing daily hourly data
;   - Computing daily averages using simple loop plus array syntax
;   - Dealing with (possibly) packed data on netCDF file
;***************************************************************
;
; These files are loaded by default in NCL V6.2.0 and newer
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"  
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"  
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"  

begin
 diri   = "./"                                 ; input dir
 fili   = systemfunc("cd "+diri+" ; ls *as_CRCM_*nc") ; input file name(s)
 nfil   = dimsizes(fili)
 print(fili)

;varName= (/"uas", "vas"/)                 ; manually specify
 varName= str_get_field(fili, 1, "_")      ; v5.1.1

 netCDF = True
 PLOT   = True

 if (netCDF) then
     ncDir   = "./"                        ; directory for netCDF output
     ncPack  = True                        ; True=>create "short"; False=>float
     if (ncPack) then
         scale   = 1000.
         optPack = True
         optPack@scale_factor = 1./scale   ; COARDS/CF convention attributes
         optPack@add_offset   = 0.
     end if
 end if

 if (PLOT) then
     pltDir  = "./"                        ; directory for plot output
     pltName = "crcm"                      ; plot name output
     pltType = "png"                       ; send graphics to PNG file  
     wks     = gsn_open_wks(pltType, pltDir+pltName)
 end if

;*******************************************
; Loop over files                                    
;*******************************************
 do nf=0,nfil-1
    f       = addfile (diri+fili(nf), "r")
    x3      = f->$varName(nf)$       ; (time,yc, xc) 3 hrly
    printVarSummary(x3)

    dimx3   = dimsizes(x3)
    NTIM    = dimx3(0)               ; TOTAL number of time steps
    print("NTIM="+NTIM)

;*******************************************
; Create daily averages                              
;*******************************************
    ntJump  = 24                     ; hourly

    x       = x3(::ntJump,:,:)       ; trick: create array with meta
    printVarSummary(x)               ; values will be overwritten with averages

    ntStrt  = 0
    ntLast  = ntJump-1
    do nt=0,NTIM-1,ntJump            ; dim_avg_n v5.1.1
       x(nt/ntJump,:,:) = (/ dim_avg_n(x3(ntStrt:ntLast,:,:), 0) /)   ; (/.../) ignore meta
       ntStrt  = ntStrt+ntJump
       ntLast  = ntLast+ntJump
    end do
    x@info_tag     = "daily average"

;*******************************************
; Miscellaneous for netCDF or graphics
;*******************************************
 
    time    = x&time                 ; time associated with daily averages
    lat     = f->lat                 ; lat(yc, xc)
    lon     = f->lon                 ; lon(yc, xc)

    filSfx  = get_file_suffix(fili(nf),0)  
    filRoot = filSfx@fBase           ; naming convenience
    print("filRoot="+filRoot)

    ymd     = cd_calendar(time,-2)   ; yyyymmdd (integer)
    yfrac   = cd_calendar(time, 4)   ; year.fraction_of_year

    ntim    = dimsizes(time)              
    print("ntim="+ntim)
       
    yyyy    = ymd/10000
    mmdd    = ymd - (yyyy*10000)
    mm      = mmdd/100
    dd      = mmdd-(mm*100)
    hh      = 12             ; center of 'mass' for the day

    yrStrt  = yyyy(0)
    yrLast  = yyyy(ntim-1)

;********************************************
; create plot
;********************************************
    if (PLOT) then
        pltNamV = pltName+"."+yrStrt+"-"+yrLast+"."+varName(nf)
;        wks     = gsn_open_wks(pltType, pltDir+pltName)
                                             
        x@lat2d = lat                 
        x@lon2d = lon
   
        res                       = True               ; plot mods desired
        res@gsnMaximize           = True               ; make ps/pdf large
        res@tiMainString          = fili(nf)
   
        res@cnFillOn              = True               ; color fill  
        res@cnFillPalette         = "amwg"             ; set color map
        res@cnLinesOn             = False              ; no contour lines
        res@cnLineLabelsOn        = False              ; no contour labels
        res@cnInfoLabelOn         = False              ; no contour info label
       ;res@cnFillMode            = "RasterFill"       ; Raster Mode

       ;res@cnLevelSelectionMode = "ManualLevels"     ; set manual contour levels
       ;res@cnMinLevelValF       = -16.               ; set min contour level
       ;res@cnMaxLevelValF       =  16.               ; set max contour level
       ;res@cnLevelSpacingF      =   2.               ; set contour spacing

        res@gsnAddCyclic          = False  
        res@mpFillOn              = False              ; no color fill
        res@mpMinLatF             = min(lat)           ; Entire grid
        res@mpMaxLatF             = max(lat)
        res@mpMinLonF             = min(lon)
        res@mpMaxLonF             = max(lon)
        
        res@gsnCenterString       = ymd(0)

        plot = gsn_csm_contour_map(wks,x(0,:,:),res)  ; 1st time only     
   end if   ; PLOT

;************************************************
; Create netCDF: Use NCL's Simple Method
;************************************************
  
   if (netCDF) then
       if (isatt(x,"lat2d")) then
           delete(x@lat2d)
           delete(x@lon2d)
       end if

       if (ncPack) then
           x = round(x*scale, 0)/scale     ; avoid truncation bias
       end if

       xMin  = min(x)
       xMax  = max(x)
       x@actual_range = (/ xMin, xMax /)  

       date     = yyyy*1000000 + mm*10000 + dd*100 
       date!0   = "time"
       date@units = "yyyymmdd"
       
       globeAtt              = 1
       globeAtt@title        = "CRCM: daily averages"
      ;globeAtt@data_source  = "..."
       globeAtt@acronym      = "CRCM: Canadian Regional Climate Model"
       globeAtt@description  = "http://www.cccma.ec.gc.ca/models/crcm.shtml"
       globeAtt@creation_date= systemfunc ("date" )

       ncVarName = str_upper(varName(nf))  ; make nc name upper case

       ncFil     = "DailyAvg."+yrStrt+"-"+yrLast+"."+ncVarName+".nc" 
     
       NCFILE    = ncDir + ncFil
       system ("/bin/rm -f " +  NCFILE)    ; remove any pre-exist file 
           
       ncdf   = addfile(NCFILE,"c")     
     
       fileattdef( ncdf, globeAtt )        ; create the global [file] attributes
       filedimdef(ncdf,"time",-1,True)     ; force "time" to be unlimited
                                         
       ncdf->time   = time
       ncdf->lat    = lat
       ncdf->lon    = lon
       ncdf->date   = date

       if (ncPack) then
           x_short = pack_values(x, "short", optPack)
           ncdf->$ncVarName$ = x_short
       else
           ncdf->$ncVarName$ = x 
       end if
   end if   ; netCDF

 end do     ; nf

end