load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"

begin

  monthdef = (/"","Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep", \
                    "Oct","Nov","Dec"/)

  analperiod  = (/"2016-09-19_06:00:00", "2017-01-27_06:00:00"/)
  level = 700   ; hPa

;------------
  DATADir = "/data/WRF/seasonal/sep20NDJ16/WRFRUN/"
  FILES = systemfunc (" ls -1 " + DATADir + "wrfout_d01*")
  a = addfiles(FILES+".nc","r")
 print(a)
;---- find index

  times = wrf_user_getvar(a,"times",-1)  ; get all times in the file
  print(times)
  ntimes = dimsizes(times)         ; number of times in the file

  x = a[:]->Times
 printVarSummary(x)
  time0 = wrf_times_c(x, 0)
   printVarSummary(time0)
  fdate = cd_calendar(time0, 0)
  fmonth  = tointeger(fdate(:,1))
  fday    = tointeger(fdate(:,2))
  timelabels = monthdef(fmonth)+"-"+fday
  delete(x)
  ihrs      = get1Dindex(times,analperiod)
  lendays    = (ihrs(1)-ihrs(0))/4
  timedays  = times(ihrs(0):ihrs(1):4)
  startindx = get1Dindex(times,timedays(0))
  endindx = get1Dindex(times,timedays(dimsizes(timedays)-1))
  print(times(startindx)+ " " + startindx + " " +times(startindx+20))  

;---Read, vertical velocity(w), u and v winds, temp(tc), pressure(p), relative humidity(RH),one by one to reduce memory use

  ;p0 = wrf_user_getvar(a,"pressure", 20  )   ; pressure
  p0 = wrf_user_getvar(a,"pressure", 20  )   ; pressure
  dims = dimsizes(p0)
  print(dims)
  ;print(zew)
  varT  = new( (/40,dims(1),dims(2)/), typeof(p0))
  ;arT = p0(0:40,1,:,:)
  ;arT = -9999
  ;varT@_FillValue = -9999
  varU = varT
  varwa = varT
  varRh = varT
  varV  = varT
  delete(p0)

do itimes = startindx, startindx+40
  p = wrf_user_getvar(a,"pressure",itimes)   ; pressure
  tc  = wrf_user_getvar(a,"tc",itimes)    ; 3D temperature
  ;tc = dim_avg_n_Wrap(tc1(startindx:endindx,:,:,:),0)
  ;delete(tc1)
  tc_plane = wrf_user_intrp3d(tc,p,"h",level,0.,False)
  delete(tc)
   printVarSummary(tc_plane)
   
  rh = wrf_user_getvar(a,"rh",itimes)        ; relative humidity
 ; rh  = dim_avg_n_Wrap(rh1(startindx:endindx,:,:,:),0)
  ;delete(rh1)
  rh_plane = wrf_user_intrp3d(rh,p,"h",level,0.,False)
  ;delete(rh)
  
  u   = wrf_user_getvar(a,"ua",itimes)    ; 3D U at mass points
  ;u  = dim_avg_n_Wrap(u1(startindx:endindx,:,:,:),0)
  ;delete(u1)
  u_plane  = wrf_user_intrp3d( u,p,"h",level,0.,False)
  ;delete(u)

  v   = wrf_user_getvar(a,"va",itimes)    ; 3D V at mass points
  ;v  = dim_avg_n_Wrap(v1(startindx:endindx,:,:,:),0)
  ;delete(v1)
  v_plane  = wrf_user_intrp3d( v,p,"h",level,0.,False)
  ;delete(v)

 wa   = wrf_user_getvar(a,"wa",itimes)    ; 3D W vertical velocity at mass points
  ;w = dim_avg_n_Wrap(w1(startindx:endindx,:,:,:),0)
  ;delete(w1)
  wa_plane  = wrf_user_intrp3d( wa,p,"h",level,0.,False)
  ;delete(w)
  printVarSummary(wa)

  varT(itimes,:,:) =  tc_plane 
  varU(itimes,:,:) =   u_plane 
  varV(itimes,:,:) =   v_plane 
  varRh(itimes,:,:) =   rh_plane 

  varwa(itimes,:,:) =   wa_plane

end do
  XLAT      = wrf_user_getvar(a,"XLAT",itimes)    ; XLAT
  XLONG      = wrf_user_getvar(a,"XLONG",itimes)    ; XLONG

 
;-----------------------saving output file-----------------------------------------
    system("/bin/rm -f  file850wind.nc")

     ncdf = addfile("file850wind.nc" ,"c")  ; open output netCDF file
;---------------------assigning units-----------------------
;      tc_plane@description = "Temperature at " + level 
;      tc_plane@units       = "degC"
;      w_plane@description = "vertical velocity at " + level
;      w_plane@units       = "m/s"
;      u_plane@units = "m/s"
;      v_plane@units = "m/s"
;      rh_plane@description = "Relative Humidity at " + level 
;      rh_plane@units       = "%"


    ;---------------------global attributes-----------------------

       fAtt               = True            ; assign file attributes

       fAtt@title         = "u-v-w winds "

       fAtt@Conventions   = "None"

       fAtt@creation_date = systemfunc ("date")

       fileattdef( ncdf, fAtt )            ; copy file attributes

        ncdf->U  = varU
        ncdf->V  = varV
        ncdf->wa  = varwa
        ncdf->T  = varT
        ncdf->RH = varRh
        ncdf->XLAT=XLAT
        ncdf->XLONG=XLONG
       
       
end