; ncl_talk: 3 Aug 2017
; We are trying to find a script that would read 24 individual 1-hour wrfout_d01* 
; output files for every day of an entire month of the year, then average the 
; individual hourly ground level ozone (o3) concentrations data over each daily 
; 24-hour period (at lowest level nl=0), and write a new netcdf file with an gridded 
; o3 average concentration for that specific day. The end result would be a visual 
; representation of a gridded o3 average concentration for any specific day, if possible. 

load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"

   dirw = "./"
;  dirw = "./wrfout_d01_2010-07-15_16:00:00.nc"
;  dirw = "./wrfout_d01_2010-07-15_17:00:00.nc"
;  dirw = "./wrfout_d01_2010-07-15_18:00:00.nc"
;  dirw = "./wrfout_d01_2010-07-15_19:00:00.nc"
;  dirw = "./wrfout_d01_2010-07-15_20:00:00.nc"
  
  filw = systemfunc("cd "+dirw+" ; ls wrfout_d01_2010-07-15*")
;  filw = systemfunc("cd "+dirw+" ; ls wrfout_d01_2010-07-16*")
;  filw = systemfunc("cd "+dirw+" ; ls wrfout_d01_2010-07-17*")
;  filw = systemfunc("cd "+dirw+" ; ls wrfout_d01_2010-07-18*")
;  filw = systemfunc("cd "+dirw+" ; ls wrfout_d01_2010-07-19*")
;  filw = systemfunc("cd "+dirw+" ; ls wrfout_d01_2010-07-15_20:00:00")
;  filw = systemfunc("cd "+dirw+" ; ls dailyave_2010-07-15.nc")
  
   pthw = dirw+filw
   print(pthw)
  
;---https://www.ncl.ucar.edu/Document/Functions/Built-in/addfiles.shtml
  f    = addfiles(pthw, "r")                    ; 'f' is type list

  Times= f[:]->Times                            ; char Times(Time, DateStrLen) 

  Times_str= tostring( Times )          
  print(Times_str)
  print("-----")

;---https://www.ncl.ucar.edu/Document/Functions/WRF_contributed/wrf_times_c.shtml

  Time = wrf_times_c( f[:]->Times,0 )          ; CF/COARDS 'time' variable
  if (isatt(Time,"_FillValue")) then
      delete(Time@_FillValue)
  end if
  print(Time)
  print("-----")

; exit 
  
;---Read 'o3' across all files

  o3   = wrf_user_getvar(f,"o3",-1)             ; o3(Time, bottom_top, south_north, west_east)
  printVarSummary(o3)
  print("-----")
  
;---Associate 'Time' coordinate variable with variables

  Times&Time = Time
  o3&Time    = Time 
  printVarSummary(Times)
  print("-----")
  printVarSummary(o3)
  print("-----")

;---https://www.ncl.ucar.edu/Document/Functions/Contributed/calculate_daily_values.shtml
;---requires 6.4.0

  klv     = 0                                   ; lowest level
  nDim    = 0         
  o3DayAvg= calculate_daily_values(o3(:,klv,:,:),"max", nDim, False) ; level 'klv' only
  
;  o3DayAvg@description = "Surface Ozone Concentration: " + o3@description
  o3DayAvg@level       = "klv="+klv+"; lowest level"
  o3DayAvg=o3DayAvg*1000
  o3DayAvg@units="ppbv"

  printVarSummary(o3DayAvg)
  printMinMax(o3DayAvg,0)
  print("Daily Average Value")
  print(avg(o3DayAvg))
  print("-----")

;---netCDF"

  diro = "./"                                   ; output directory
  filo = "o3.day.nc"
  ptho = diro + filo
  system ("/bin/rm -f "+ptho)
  fnc = addfile (ptho, "c")
  
  filAtt = 0
  filAtt@title = "O3: Daily Max"
  filAtt@Conventions = "None"
  filAtt@creation_date = systemfunc("date")
  filAtt@NCL    = get_ncl_version()
  
  copy_VarAtts(f[0], filAtt)                   ; copy WRF file attributes 
 
  fileattdef(fnc, filAtt )                     ; create file attributes
  filedimdef(fnc,"Time",-1,True)               ; make Time UNLIMITED
  
;  fnc->Time     = Time 
;  fnc->Times    = Times                       ; original character
;  fnc->XLAT     = f[0]->XLAT(0,:,:)           ; no 'Time' dimension
;  fnc->XLONG    = f[0]->XLONG(0,:,:)
;  fnc->o3DayAvg = o3DayAvg 
  
  fnc->o3DayAvg = o3DayAvg                     ; write meta data associated with daily avg
  fnc->XLAT     = f[0]->XLAT(0,:,:)            ; no 'Time' dimension
  fnc->XLONG    = f[0]->XLONG(0,:,:)