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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "./shapefile_mask_data.ncl"

begin


;--- Years of interest
  YS = 1990 						; Which years are we interested in
  YE = 2000
  numYEARS = YE-YS+1					; How many years needs to be processed

  YEARS = new (numYEARS, integer)			; Create an array containing years of interest
  YEARS!0 = "Year"
  do iyy = 0,numYEARS-1					
    YEARS(iyy) = YS + iyy
  end do

;--- Input Data
  ensembles = (/"rtty_1990"/)		; Name of the regional climate model ensemble we are using
  numENS = dimsizes(ensembles)				; How many ensembles do we have
  ;DIR = "/kumquat/rcr/rctutorial_2016/Day1/data/"	; In which directory is the input data
  DIR = ""
  rootFILE = "t2_maxmin_daily_jja_"			; Root name of the input files we are going to read

  wrf_filename  = DIR+rootFILE+YEARS(0)+"_"+ensembles(0)+".nc"		; open files
  a = addfile(wrf_filename,"r")


;--- Shapefiles and Areas we are interested in
  ;DIRs = "../data/ShapeFiles_US/"					; Location of shapefile
  DIRs = "data/ShapeFiles_US/"					; Location of shapefile
  shp_filename = DIRs + "cb_2014_us_state_20m.shp"		; Shapefile names

  regions = (/"Colorado"/) ; Our areas of interest. Need to know which areas are available in the Shapefiles. 
  numRegions = dimsizes(regions)			; How many Regions are we going to process
  DataNames = new ( (/2,numRegions/) , string)		; Create an array where we are going to store some new variable names

; Loop over the areas we are interested in and create the new variables names (i.e., T2 Min and Max for each area)
  do ireg = 0,numRegions-1
    foo1 = stringtochar(shp_filename)
    ind1 = str_index_of_substr(shp_filename , "/", -1) + 1
    ind2 = str_index_of_substr(shp_filename , ".", -1) - 1
    foo2 = foo1(ind1:ind2)
    DataNames(0,ireg) = "T2MIN_"+regions(ireg)
    DataNames(1,ireg) = "T2MAX_"+regions(ireg)
    delete( [/foo1,foo2/] )
  end do


;--- Loop over all ensembles of interest and extract the data over the area of interest
  do iens = 0,numENS-1
    print("Working on ensemble " + ensembles(iens) )
    outfile = "t2_daily_Regional_Data_US_jja_"+ensembles(iens)+".nc"	; Create a new output file
    out = addfile(outfile,"c")
    dim_names = (/ "Year", "JulianDay" /)
    dim_sizes = (/  -1, 365 /)
    dimUnlim  = (/  True, False /)
    filedimdef( out, dim_names, dim_sizes, dimUnlim )			; The output file will have dimensions of Years (open-ended dimension) and Days (365 or 93) per season
    out->YEARS = YEARS							; Write years information to the new output file


    do iyy = 0,numYEARS-1						; Loop over all years
      print("   Working on year " + YEARS(iyy) )
      wrf_filename  = DIR+rootFILE+YEARS(iyy)+"_"+ensembles(iens)+".nc"	; Open files
      a = addfile(wrf_filename,"r")

      T2MIN       = a->T2MIN				; Read T2 MIN and MAX from the file and make sure we also get the associated lat and lon attached to these variables
      T2MIN@lat2d = a->XLAT
      T2MIN@lon2d = a->XLONG
      T2MAX       = a->T2MAX
      T2MAX@lat2d = a->XLAT
      T2MAX@lon2d = a->XLONG
      dims = dimsizes(T2MAX)				; How many times are in the file (i.e., days - we know it is 93, but best to not hardcode)
      numDays = dims(0)

  
      opt             = True				; Set options to read Shapefile
      opt@shape_var   = "NAME"				; We know the shape names are stored in the array NAME inside the shapefiles
      opt@debug       = False

      do ireg = 0,numRegions-1				; Loop over regions
        print("      Region " + regions(ireg) ) 
        opt@shape_names = regions(ireg)			; Make sure to specific the region we are working with

        do iday = 0,numDays-1				; Loop over days of the year
          
          var_mask  = shapefile_mask_data(T2MIN(iday,:,:),shp_filename,opt) 	; Read T2MIN and set all data values to missing except for those over our region of interest
          var_mask_1d = ndtooned(var_mask)					; Turn this masked dataset into a 1d array

          if (iday .eq. 0) then							; Create new arrays - only need to do this once per area of interst
            varRegIDs = ind(.not.ismissing(var_mask_1d))
            varReg = var_mask_1d(varRegIDs)
            dimReg = dimsizes(varReg)
            delete(varReg)

            min_region = new ( (/365,dimReg/), float)
            min_region!0 = "JulianDay"
            min_region!1 = "datadim"+ireg
            min_region@description = "DAILY MIN TEMP at 2 M for the region"
            min_region@units = "C"

            max_region = new ( (/365,dimReg/), float)
            max_region!0 = "JulianDay"
            max_region!1 = "datadim"+ireg
            max_region@description = "DAILY MAX TEMP at 2 M for the region"
            max_region@units = "C"

            if (iyy .eq. 0) then						; Add new dimensions to the output data file
              dd = "datadim"+ireg
              filedimdef( out, dd, dimReg, False )
              DimNames = (/ "Year", "JulianDay", dd /)
              new_var_type = new ( 2 , string )
              new_var_type = "float"
              filevardef(out, DataNames(:,ireg), new_var_type, DimNames )
            end if
          end if

          min_region(iday,:) = (/var_mask_1d(varRegIDs)/)				; Save the data to our new array
          delete([/var_mask,var_mask_1d/])

          var_mask  = shapefile_mask_data(T2MAX(iday,:,:),shp_filename,opt) 		; Do the same for T2 Max
          var_mask_1d = ndtooned(var_mask)			
          max_region(iday,:) = (/var_mask_1d(varRegIDs)/)
          delete([/var_mask,var_mask_1d/])

        end do		; End loop of days


        ;--- Write out the data to the file. Note that we do this per year, and we can because year does not have a fixed dimension.
        out->$DataNames(0,ireg)$(iyy,:,:) = min_region
        out->$DataNames(1,ireg)$(iyy,:,:) = max_region
        delete([/varRegIDs,min_region,max_region/])
      
      end do		; End of region loop


      delete ( [/T2MIN, T2MAX/] )

    end do		; End of year loop

  end do		; End of ensemble loop

end