; 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"
load "./calculate_monday_values.ncl"
   dirw = "./"

;  filw = systemfunc("cd "+dirw+" ; ls wrfout_d01_2010-01*")
;  filw = systemfunc("cd "+dirw+" ; ls wrfout_d01_2010-04*")
;  filw = systemfunc("cd "+dirw+" ; ls wrfout_d01_2010-05*")
;  filw = systemfunc("cd "+dirw+" ; ls wrfout_d01_2010-06*")
   filw = systemfunc("cd "+dirw+" ; ls wrfout_d01_2010-07*")
;  filw = systemfunc("cd "+dirw+" ; ls wrfout_d01_2010-[07..08]*")
;  filw = systemfunc("cd "+dirw+" ; ls wrfout_d01_2010-08*")
;  filw = systemfunc("cd "+dirw+" ; ls wrfout_d01_2010-09*")
;  filw = systemfunc("cd "+dirw+" ; ls wrfout_d01_2010-10*")

  
   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         
  o3MonthAvg= calculate_monthly_values(o3(:,klv,:,:),"avg", nDim, False) ; level 'klv' only
  
;  o3MonthAvg@description = "Surface Ozone Concentration: " + o3@description
  o3MonthAvg@level       = "klv="+klv+"; lowest level"
  o3MonthAvg=o3MonthAvg*1000
  o3MonthAvg@units="ppbv"

  printVarSummary(o3MonthAvg)
  printMinMax(o3MonthAvg,0)
  print("Monthly Average Value")
  print(avg(o3MonthAvg))
  print("-----")

;---netCDF"

  diro = "./"                                   ; output directory
  filo = "o3.month.nc"
  ptho = diro + filo
  system ("/bin/rm -f "+ptho)
  fnc = addfile (ptho, "c")
  
  filAtt = 0
  filAtt@title = "O3: Monthly Avg"
  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 = o3MonthAvg 
  
  fnc->o3MonthAvg = o3MonthAvg                     ; write meta data associated with daily avg
  fnc->XLAT     = f[0]->XLAT(0,:,:)            ; no 'Time' dimension
  fnc->XLONG    = f[0]->XLONG(0,:,:)