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

; Ensemble member(s) to be processed
mem = (/"rtty_1990"/)
numMEM = dimsizes(mem)

; Years we are interested in
years = (/"1990","2003"/)
numYEARS = dimsizes(years)

; Loop through the number of ensemble members
do k=0,numMEM-1

; Loop through the number of years
  do j=0,numYEARS-1

    print("Year = " +years(j))
; Create the netdcf output file
    out = addfile("t2_maxmin_daily_jja_" + years(j) + "_" + mem(k) + ".nc","c")
    filedimdef(out,"Time",-1,True)

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

; Create new variables for the daily max and min for both masked and unmasked
    tmaxday = new((/numFILES,279,395/),float)
    tlandmaxday = new((/numFILES,279,395/),float)
    tminday = new((/numFILES,279,395/),float)
    tlandminday = new((/numFILES,279,395/),float)

; Because wrf is sometime weird, there is a problem with the landmask in the output
;   files.  Therefore we are taking the landmask from a met_em wrf input file.
    a = addfile("data/met_em.d01.2004-12-16_12:00:00.nc","r")
    landmask = a->LANDMASK(0,:,:)

; Looping through a years 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) then
        fileattdef(out,a)
      end if

; Read in temperature, latitude and longitude
      t = a->T2-273.15    ; We want t2 in Celcius!
      lat = a->XLAT(0,:,:)
      lon = a->XLONG(0,:,:)

; Masking out the oceans using the landmask to just get t2 over land
      tland = mask(t,landmask,1)

; Calculated maxes for t2 and t2mask
      tmax = dim_max_n(t,0)
      tlandmax = dim_max_n(tland,0)

; Writing the maxes to the maxday variables
      tmaxday(i,:,:) = tmax
      tlandmaxday(i,:,:) = tlandmax

; Doing the same for minimums...
      tmin = dim_min_n(t,0)
      tlandmin = dim_min_n(tland,0)

      tminday(i,:,:) = tmin
      tlandminday(i,:,:) = tlandmin

; Deleting t and tland to start fresh!
      delete([/t,tland/])

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

; Writing out lat, lon, and landmask to the netcdf file
    out->XLAT = lat
    out->XLONG = lon
    out->LANDMASK = landmask

; Writing out the variable attributes and values to the netcdf output file
    tmaxday!0 = "Time"
    tmaxday!1 = "south_north"
    tmaxday!2 = "west_east"
    tmaxday@units = "C"
    tmaxday@coordinates = "XLONG XLAT"
    tmaxday@description = "DAILY MAX TEMP at 2 M"
    out->T2MAX = tmaxday

    tlandmaxday!0 = "Time"
    tlandmaxday!1 = "south_north"
    tlandmaxday!2 = "west_east"
    tlandmaxday@units = "C"
    tlandmaxday@coordinates = "XLONG XLAT"
    tlandmaxday@description = "DAILY MAX TEMP at 2 M (masked)"
    out->T2MAX_MASK = tlandmaxday

    tminday!0 = "Time"
    tminday!1 = "south_north"
    tminday!2 = "west_east"
    tminday@units = "C"
    tminday@coordinates = "XLONG XLAT"
    tminday@description = "DAILY MIN TEMP at 2 M"
    out->T2MIN = tminday

    tlandminday!0 = "Time"
    tlandminday!1 = "south_north"
    tlandminday!2 = "west_east"
    tlandminday@units = "C"
    tlandminday@coordinates = "XLONG XLAT"
    tlandminday@description = "DAILY MIN TEMP at 2 M (masked)"
    out->T2MIN_MASK = tlandminday

    delete([/FILES,numFILES,a,out,tmax,tmaxday,tlandmaxday,tmin,tminday,tlandminday/])

; End of the loop over all the years we care about
  end do

; End of the loop over ensemble members
end do

; End of program
end