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/wrf/WRFUserARW.ncl"

begin

; Years and months we are interested in
years = (/"1990","1991","1992","1993","1994","1995","1996","1997","1998","1999","2000"/)
numYEARS = dimsizes(years)

mons = (/"06","07","08"/)
numMONS = dimsizes(mons)

; Make new variables for diurnal calculations from the model
phour = new((/24,279,395/),float)
pocchour = new((/24,279,395/),float)
pinthour = new((/24,279,395/),float)
phouryr = new((/numYEARS,24,279,395/),float)
pocchryr = new((/numYEARS,24,279,395/),float)

; Create the netdcf output file
  out = addfile("prec_diurnal_jja_rtty.nc","c")
  filedimdef(out,"Time",-1,True)

; Loop through the number of years
  do j=0,numYEARS-1
    phourmon = new((/numMONS,24,279,395/),float)
    pocchrmon = new((/numMONS,24,279,395/),float)
    print(years(j))

; Loop through the months
    do l=0,numMONS-1

; Find the data files for the year and month of interest, list them,
;   and put them in an array of filenames to be read
    DATADir = "/sysdisk1/C3WE/tutorial/model_data/"
    FILES = systemfunc (" ls -1 " + DATADir + "auxhist1_d01_" + years(j) + "-" + mons(l) + "* ")
    numFILES = dimsizes(FILES)

; Create new variables for hourly precip amount and occurence for each hour
    phourday = new((/24,numFILES,279,395/),float)
    pocchrday = new((/24,numFILES,279,395/),float)
    pocchrday = 0

; Looping through a months worth of files
    do i=0,numFILES-1

; Reading in a day at a time and moving the file attributes to our output file
      a = addfile(FILES(i)+".nc","r")
      if(i.eq.0 .and. l.eq.0 .and. j.eq.0) then
        fileattdef(out,a)
      end if

; Here we are pulling out the specific year, month, and day from each file name.
      year = str_get_cols(FILES(i),48,51)
      mon = str_get_cols(FILES(i),53,54)
      day = str_get_cols(FILES(i),56,57)
      ;year = str_get_cols(FILES(i),24,27)
      ;mon = str_get_cols(FILES(i),29,30)
      ;day = str_get_cols(FILES(i),32,33)

; Because our precip is accumulated we need to set the next day for the end of each month.
      if(mon.eq."06" .and. day.eq."30") then
        NEWFILE = systemfunc (" ls -1 " + DATADir + "auxhist1_d01_"+year+"-07-01_00*")
        b = addfile(NEWFILE+".nc","r")
      else if(mon.eq."07" .and. day.eq."31") then
        NEWFILE = systemfunc (" ls -1 " + DATADir + "auxhist1_d01_"+year+"-08-01_00*")
        b = addfile(NEWFILE+".nc","r")
      else if(mon.eq."08" .and. day.eq."31") then
        NEWFILE = systemfunc (" ls -1 " + DATADir + "auxhist1_d01_"+year+"-09-01_00*")
        b = addfile(NEWFILE+".nc","r")
      else
        b = addfile(FILES(i+1)+".nc","r")
      end if
      end if
      end if

; Read in precipitation, latitude and longitude

      lat = a->XLAT
      lon = a->XLONG

; Loop through each hour to calculate hourly precipitation.
      do m=0,23
        if(m.eq.23) then
          rainc = a->RAINC(m,:,:)
          rainnc = a->RAINNC(m,:,:)
          irainc = a->I_RAINC(m,:,:)    ; irainc, irainnc, etc are bucket accumulations
          irainnc = a->I_RAINNC(m,:,:)
          rainc2 = b->RAINC(0,:,:)
          rainnc2 = b->RAINNC(0,:,:)
          irainc2 = b->I_RAINC(0,:,:)
          irainnc2 = b->I_RAINNC(0,:,:)
        else
          rainc = a->RAINC(m,:,:)
          rainnc = a->RAINNC(m,:,:)
          irainc = a->I_RAINC(m,:,:)    ; irainc, irainnc, etc are bucket accumulations
          irainnc = a->I_RAINNC(m,:,:)
          rainc2 = a->RAINC(m+1,:,:)
          rainnc2 = a->RAINNC(m+1,:,:)
          irainc2 = a->I_RAINC(m+1,:,:)
          irainnc2 = a->I_RAINNC(m+1,:,:)
        end if
        precc = (irainc2-irainc)*100 + (rainc2-rainc)
        precnc = (irainnc2-irainnc)*100 + (rainnc2-rainnc)
; Add convective and non-convective precipitation to get total
        prec = precc+precnc
        phourday(m,i,:,:) = prec
; Find if an occurence, if not set precip to zero
        pocchrday(m,i,:,:) = where(phourday(m,i,:,:).ge.2.5,pocchrday(m,i,:,:)+1,pocchrday(m,i,:,:))
        phourday(m,i,:,:) = where(phourday(m,i,:,:).ge.2.5,phourday(m,i,:,:),0)
        delete([/rainc,irainc,rainnc,irainnc,precc,precnc,prec/])
      end do

; End of the loop over the months worth of files
    end do

; Average the hourly precip over each month.
    phourmon(l,:,:,:) = dim_sum_n(phourday,1)
    pocchrmon(l,:,:,:) = dim_sum_n(pocchrday,1)

    delete([/FILES,numFILES,NEWFILE,b,a,phourday,pocchrday/])

; End of the monthly loop
  end do

; Average the hourly precip over each year
  phouryr(j,:,:,:) = dim_sum_n(phourmon,0)
  pocchryr(j,:,:,:) = dim_sum_n(pocchrmon,0)

  delete([/phourmon,pocchrmon/])
;End of the loop over all the years we care about
end do

; Average the precip so we just have one diurnal cycle and calculate the intensity.
pocchryr = where(pocchryr.eq.0,0.00001,pocchryr)
phour = dim_avg_n(phouryr,0)
pocchour = dim_avg_n(pocchryr,0)
phrtest = where(phour.eq.0,phour@_FillValue,phour)
pinthour = phrtest/pocchour

; Write out data to new netcdf file.
out->XLAT = lat
out->XLONG = lon
phour!0 = "Time"
phour!1 = "south_north"
phour!2 = "west_east"
phour@units = "mm"
phour@coordinates = "XLONG XLAT"
phour@description = "HOURLY TOTAL PRECIPITATION OVER JJA"
out->PREC = phour
pocchour!0 = "Time"
pocchour!1 = "south_north"
pocchour!2 = "west_east"
pocchour@coordinates = "XLONG XLAT"
pocchour@description = "NUMBER OF OCCURENCES ABOVE 2.5MM OVER JJA"
out->PRECOCC = pocchour
pinthour!0 = "Time"
pinthour!1 = "south_north"
pinthour!2 = "west_east"
pinthour@units = "mm/hour"
pinthour@coordinates = "XLONG XLAT"
pinthour@description = "INTENSITY OF OCCURENCES ABOVE 2.5MM OVER JJA"
out->PRECINT = pinthour

; End of program
end