;;<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<;;

load "$NCARG_ROOT/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/nclscripts/csm/shea_util.ncl"
load "/home/op/ventrice/nclscripts/home/ut_string.ncl"
load "/home/op/ventrice/nclscripts/ct13/scripts/comp_xy.ncl"
load "/home/op/ventrice/nclscripts/ct13/scripts/shadecomp_xy.ncl"
load "/home/op/ventrice/nclscripts/ct13/scripts/sigshade_xy.ncl"
load "/home/op/ventrice/nclscripts/ct13/scripts/signif_xy.ncl"
load "/home/op/ventrice/nclscripts/ct13/scripts/comp_xy.ncl"
load "/home/op/ventrice/nclscripts/home/tTick.ncl"
;;<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<;;

;;\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\;;
;;>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>;;
;;>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>;;
;;///////////////////////////////////////////////////////////////////////////;;
begin

  system( "date" )
  print( "Here we go!" )


 year = stringtointeger( systemfunc( "date +'%Y' -d today " ) )
  month = stringtointeger( systemfunc( "date +'%m' -d today " ) )
  day = stringtointeger( systemfunc( "date +'%d' -d today " ) )
  hour = stringtointeger( systemfunc( "date +'%H' -d today " ) )
  inHH = hour

 
timeUnits = "hours since 1800-01-01 00:00:00"

 
 if( ( inHH.ge.12 ).and.( inHH.lt.18 ) ) then
    HH = "06"
 end if

if( ( inHH.ge.18 ).or.( inHH.lt.0 ) ) then
    HH = "12"
end if

if( ( inHH.ge.0 ).and.( inHH.lt.6 ) ) then
    HH = "18"
  DD_cov = stringtoint(DD)
  DD_new = DD_cov - 1

end if

if( ( inHH.ge.6 ).and.( inHH.lt.12 ) ) then
    HH = "00"
end if


   currDate = ut_inv_calendar( year, month, day, hour, 0, 0, timeUnits, 0 )
    lastDate = ut_inv_calendar( year, month, day, hour, 0, 0, timeUnits, 0 ) + 15
    


  print(ut_string(currDate, ""))
  print(ut_string(lastDate, ""))

;---------------------------------------------------------------


 pltTitle = month+"/"+day+"/"+year+"/"+inHH+"Z "  

  plotType = "png"
  plotDir = "/home/op/ventrice/real_time/figures/"
  plotName = plotDir+"Forecast_ECM_USA_SUMMER"

   fontHeightF = 0.012


 ;bin1
  latRange =  (/35, 50 /)          
  lonRange =  (/260,285 /)   

 lon1 = 360 - lonRange(0)
 lon2 = 360 - lonRange(1)

  print( "Reading the archive mean/sigma..." )

   fileName = "/archive/data/ncep/Cropli/NaGaDex_USA_SUMMER_Sigma_Mean.nc"
  inFile = addfile( fileName, "r" )

    sigma = inFile->WiND_sigma
    mean = inFile->WiND_mean

;----------------------------

  ;Edit number of forecast days
  nfdays = 15
  model_data = nfdays
  model_data = nfdays * 24
  

  ;Time

    
  firstTime = currDate 
  lastTime = currDate + nfdays
  duration = lastTime - firstTime
  nTimes = 1 + doubletoint( duration / 6 )
  time = firstTime + fspan( 0, duration, nTimes )

  time@units = timeUnits
 
;  print(ut_string(lastTime, ""))


;------------------------ Read GFS ENS Data ------------------------



  print( "Reading Forecast Data..." )
  fHHH = ispan( 0, model_data, 6 ) 


  hr_arr = (/"00","06","12","18","24","30","36","42","48","54","60","66", \\
                "72","78","84","90","96","102","108","114", "120","126", \\
                "132", "138", "144", "150", "156", "162", "168",  "174", "180", \\
                "186", "192", "198","204", "210","216","222","228", "234","240", \\
               "246", "252", "258", "264", "270", "276", "282", "288", "294", "300", \\
                "306", "312", "318", "324", "330", "336", "342", "348", "354", "360", \\
                 "372", "378", "384" /)

    mem = (/"01","02","03","04","05","06","07","08","09","10","11","12","13","14","15","16","17","18","19","20"/)


   fileNames = new((/64, 20/), "string")
 

  do hLag = 0, dimsizes(hr_arr)-1
        do    x = 0, dimsizes(mem)-1

    fileNames(hLag,x) = "/archive/real_time/gfs/gfs_ens_fcst/"+hr_arr(hLag)+"/ens."+mem(x)

        end do
  end do

   printVarSummary(fileNames)
  ; print(fileNames(0,:))

  print("Adding files together")

   v500in = new((/64, 20,73,144/),"float")

  do y = 0,19
  gfsFiles = addfiles( fileNames(:,y) + ".grib", "r" )

  print(gfsFiles)

  ListSetType( gfsFiles, "join" )

   print("Finished joining files")


  printVarSummary(gfsFiles)


  presLev = 500 
   
  print("Constructing forecast variable")

   test = gfsFiles[:]->V_GRD_2_ISBL
   print(test(:,0,0,0))

   printVarSummary(test)

   printVarSummary(v500in(:,0,:,:))

   v500in(:,y,:,:) = gfsFiles[:]->V_GRD_2_ISBL(:,{500},::-1,:)

   printVarSummary(v500in)

end do

   printVarSummary(v500in)

exit