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

  ensembles = (/"rtty"/)
  numENS = dimsizes(ensembles)
  regions = (/"Colorado"/)
  numRegions = dimsizes(regions)

  ; Place of data files and years of interst
  DIR = ""
  rootFILE = "t2_maxmin_daily_jja_"
  rootFILE2 = "t2_maxmin_daily_jja_"
  YS = 1990 
  YE = 2000
  numYEARS = YE-YS+1
  YEARS = new (numYEARS, integer)
  YEARS!0 = "Year"
  do iyy = 0,numYEARS-1							; make an array with years
    YEARS(iyy) = YS + iyy
  end do

  wrf_filename = DIR+rootFILE+YEARS(0)+"_ck6m_CFSR.nc"      ; open files
  a = addfile(wrf_filename,"r")
  mdims = getfilevardimsizes(a,"lat_0") ; get some dimension sizes for the file
  nd = dimsizes(mdims)

  ; Place of shapefiles
  DIRs = "data/ShapeFiles_US/"
  shp_filename = DIRs + "cb_2014_us_state_20m.shp"
  opt             = True
  opt@lon360      = True
  opt@debug       = False
  opt@shape_var   = "NAME"
  DataNames = new ( (/2,numRegions/) , string)
  do ireg = 0,numRegions-1						; loop over regions and get some region information 
    opt@shape_names = regions(ireg)
    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


  do iens = 0,numENS-1							; loop over all ensembles
    print("Working on ensemble " + ensembles(iens) )
    outfile = "t2_daily_State_Data_US_"+ensembles(iens)+"_CFSR.nc"		; create new output file
    out = addfile(outfile,"c")
    dim_names = (/ "Year", "JulianDay" /)
    dim_sizes = (/  -1, 93 /)
    dimUnlim  = (/  True, False /)
    filedimdef( out, dim_names, dim_sizes, dimUnlim )
    out->YEARS = YEARS


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

      T2MIN       = a->T2MIN
      T2MIN@lat2d = a->lat_0
      T2MIN@lon2d = a->lon_0
      T2MAX       = a->T2MAX
      T2MAX@lat2d = a->lat_0
      T2MAX@lon2d = a->lon_0
      dims = dimsizes(T2MAX)
      numTimes = dims(0)
      leapyear = False							; keep track of leap years
  
      do ireg = 0,numRegions-1						; loop over regions
        opt@shape_names = regions(ireg)
        CalcArea = True
        foo  = shapefile_mask_data(T2MIN(0,:,:),shp_filename,opt)
        foo1 = ndtooned(foo)                            ; turn into 1d array
          if ( all(ismissing(foo1)) ) then
            CalcArea = False
            print("Nothing to do for " + regions(ireg) )
          end if
        delete( [/foo,foo1/] )

        if (CalcArea) then
        print("      Region " + regions(ireg) ) 
        iday = -1
        do it = 0,numTimes-1						; loop over days of the year
          
          if ( leapyear .and. it .eq. 59) then				; for leap years Feb 29 is day 60 so in NCL time skip day 59
              ;  -- We have a leap year, skipping Feb 29"
          
          else 
            iday = iday + 1 
            var_mask  = shapefile_mask_data(T2MIN(it,:,:),shp_filename,opt) ;Set all var values to missing except for those over our region of interest
            var_mask_1d = ndtooned(var_mask)				; turn into 1d array
            if (it .eq. 0) then						; Create new arrays 
              varRegIDs = ind(.not.ismissing(var_mask_1d))
              varReg = var_mask_1d(varRegIDs)
              dimReg = dimsizes(varReg)
              delete(varReg)
              min_region = new ( (/93,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 ( (/93,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
                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)/)		; pull out the data that we are intersted in and save it
            delete([/var_mask,var_mask_1d/])

            var_mask  = shapefile_mask_data(T2MAX(it,:,:),shp_filename,opt) 
            var_mask_1d = ndtooned(var_mask)			
            max_region(iday,:) = (/var_mask_1d(varRegIDs)/)
            delete([/var_mask,var_mask_1d/])

          end if	; end check for leap years
        end do	; end loop of days


        out->$DataNames(0,ireg)$(iyy,:,:) = min_region
        out->$DataNames(1,ireg)$(iyy,:,:) = max_region
        delete([/varRegIDs,min_region,max_region/])
      
      end if
      end do	; --- end of region loop

      delete ( [/T2MIN, T2MAX/] )

    end do	; --- end of year loop

  end do	; --- end of ensemble loop

end