;;>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>;;
;; make_0.5_cfsr.ncl
;; Michael Ventrice (Michael.Ventrice@weather.com)
;; Nov 2013

;;---------------------------------------------------------------------------;;
;; Requires: "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;;---------------------------------------------------------------------------;;
;; Prototype:
;;
;;         function make_raw_era_interim(
;;                  fil_in:string,
;;                  var_in:string,
;;                  var_type:integer
;;                  plev:integer
;;                  fil_out:string,
;;                  var_out:string,
;;                  long_name:string,
;;                  units:string,
;;                  nHarm:integer,
;;                  years:integer
;;         )
;;
;;---------------------------------------------------------------------------;;

;;Return Value
;;
;;A continuous NETCDF file of the specified variable is output in the specifed
;;directory
;;---------------------------------------------------------------------------;;
;;Description
;;
;;Reads in one year at a time of GRIB format era interim data and creates a
;;time continuous NETCDF file.
;;---------------------------------------------------------------------------;;

;;>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>;;
load "$NCARG_ROOT/nclscripts/csm/contributed.ncl"
;;>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>;;

undef ("make_raw_era_interim")
procedure make_raw_era_interim (dir_in:string,var_type:integer,plev:integer,fil_out:string,\
                                var_out:string,long_name:string, units:string)
begin

  months  = sprinti("%0.2i",ispan(1,12,1))
  years   = sprinti("%0.4i",ispan(1981,2010,1))

  nyears  = dimsizes(years)
  nmonths = dimsizes(months)

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

  print("Constructing a continuous NETCDF file for "+var_out+".")

;;---------------------------------------------------------------------------;;
;;Specify input files

  nfiles = nmonths*nyears

  all_fil_u = new((/nfiles/),"string")   ; Specify input file strings

  do yr = 0,nyears-1,1         
     do m = 0,nmonths-1,1
       all_fil_u(nmonths*yr+m) = dir_in+"wnd10mx0.5.gdas."+years(yr)+months(m)+".grb2"
     end do
  end do

  print("File strings created.")

;;---------------------------------------------------------------------------;;
;;Spefify metadata

  ;Lat/Lon

  lat = fspan(-90,90,361)
  lon = fspan(-180,179.5,720)

  lat@units              = "degrees_north"
  lat@long_name          = "latitude"
  lat@grid_resolution    = "0.5_degrees"
  lat@mapping            = "cylindrical_equidistant_projection_grid"
  lat@coordinate_defines = "center" 
  lat@delta_y            = 0.5
  lat@actual_range       = (/-90.,90./)

  lon@units              = "degrees_east"
  lon@long_name          = "longitude"
  lon@grid_resolution    = "0.5_degrees"
  lon@mapping            = "cylindrical_equidistant_projection_grid"
  lon@coordinate_defines = "center" 
  lon@delta_x            = 0.5
  lon@actual_range       = (/-180.,179.5/)

  ;Time

  timeUnits = "hours since 1800-01-01 00:00:00"
    
  firstTime = ut_inv_calendar( 1981, 01, 01, 00, 0, 0, timeUnits, 0 )
  lastTime = ut_inv_calendar( 2010, 12, 31, 18, 0, 0, timeUnits, 0 )
  duration = lastTime - firstTime
  nTimes = 1 + doubletoint( duration / 6 )
  time = firstTime + fspan( 0, duration, nTimes )

  time@units = timeUnits
  time@long_name = "initial time"
  time@delta_t = "0000-00-00 06:00:00" 
  time@first_time = "1981-01-01 00:00:00"
  time@last_time = "2010-12-31 18:00:00"

   forecast_time = fspan(1,6,6) ; 30* 6 or just 6?
   forecast_time@long_name = " Forecast offset from initial time"
   forecast_time@units  =  "hours"

  
;;---------------------------------------------------------------------------;;
;;Initialize the NETCDF file

  print("Setting up the output file..."+systemfunc("date"))

  ; Set file options

  setfileoption("nc","Format","NetCDF4Classic")
  setfileoption("nc","CompressionLevel",1)

  ; Specify file attributes

  fAtt               = True            ; assign file attributes

  fAtt@source_file   = dir_in+"/wnd*.grb2"
  fAtt@creation_date = systemfunc ("date")        
  fAtt@source        = "http://rda.ucar.edu/datasets/ds093.1/" 
  fAtt@references    = "Compo et al. 2010. " 
  fAtt@description   = "CFSR "+var_out
  fAtt@comments      = "Michael Ventrice (Michael.Ventrice@weather.com) " 

  ;Create file

  system("/bin/rm -f "+fil_out)  
  outFile = addfile( fil_out, "c" )

  fileattdef( outFile, fAtt )            ; Set file attributes   

  ;Specify dimension coordinates

  dimNames = (/ "time", "forecast_time","lat", "lon" /)
  dimSizes = (/ dimsizes(time), dimsizes(forecast_time),dimsizes(lat), dimsizes(lon) /)
  dimUnlim = (/ True, False, False, False /)
  
  filedimdef( outFile, dimNames, dimSizes, dimUnlim )

  filevardef( outFile, "time", typeof(time), "time" )
  filevardef( outFile, "forecast_time", typeof(forecast_time), "forecast_time" )
  filevardef( outFile, "lat", typeof(lat), "lat" )
  filevardef( outFile, "lon", typeof(lon), "lon" )
  filevardef( outFile, var_out, "float", (/ "time", "forecast_time","lat", "lon" /) )

  filevarattdef( outFile, "time", time )
  filevarattdef( outFile, "forecast_time", forecast_time )
  filevarattdef( outFile, "lat", lat )
  filevarattdef( outFile, "lon", lon )

  outFile->time = (/ time /)
  outFile->forecast_time = (/ forecast_time /)
  outFile->lat = (/ lat /)
  outFile->lon = (/ lon /)

;;---------------------------------------------------------------------------;;
;;Begin adding to the NETCDF file one month at a time

  do filn = 0,nfiles-1

    print("Grabbing the variable from "+all_fil_u(filn)+" on "+systemfunc("date"))

    u_p = addfile(all_fil_u(filn),"r")


    u = lonFlip(u_p->UGRD_P0_L103_GLL0(:,:,::-1,:))
      

   ; delete(v)
  
    ; Set metadata

    u!0 = "time"
    u!1 = "forecast_time"
    u!2 = "lat"
    u!3 = "lon"

    delete(u&lat)
    delete(u&lon)
    

    delete_VarAtts(u,-1)
  
    u@long_name  = long_name
    u@units      = units
    if (var_type.EQ.1) then   ; Pressure level or 3D cfsr variables
      u@pressure_level = plev
    end if
    u@grid_resolution = "0.5_degrees"
    u@mapping         = "cylindrical_equidistant_projection_grid"
    u@_FillValue = -999.9

    u&lat = lat
    u&lon = lon

    delete_VarAtts(u&time,-1)
    copy_VarAtts(time,u&time)


  ;print(dimsizes(time)+" "+dimsizes(forecast_time)+" "+dimsizes(lat)+" "+dimsizes(lon))
  ;printVarSummary(u)
  ;ah = dimsizes(u&time) * 30
  ;print(ah)

    outFile->$var_out$({u&time},:,:,:) = (/u/)

    if(filn.eq.0 ) then
      filevarattdef( outFile, var_out, u )  ; Not sure why this is necessary but works
    end if

    ;delete(v_p)
    delete(u_p)
    delete(u)

  end do

 



end




;;>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>;;
;;Program Main Body
;;>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>;;

begin

  print("Program begun on "+systemfunc ("date")+".")     

  dir_in  = (/"/archive/data/cfsr/1deg/wind/"/)
  plev    =  10
  var_type  = 0      ; 0 = surface, 1 = pressure level

  var_out = "u"
  fil_out = "/archive/data/cfsr/1deg/"+var_out+"."+plev+"m.nc"

  units     = "m/s"
  long_name = "10m_Uwnd"

  make_raw_era_interim(dir_in,var_type,plev,fil_out,var_out,long_name,units)

  print("Program completed on "+systemfunc ("date")+".")   

end