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  = (/"1988-10-01_06:00:00", "1988-10-15_18:00:00"/)
  level = 200.   ; hPa

;------------
  DATADir = "/home/anyuola/WRF/BML/DRY/"
  FILES = systemfunc (" ls -1 " + DATADir + "wrfout_d01_198*")
  a = addfiles(FILES+".nc","r")

;---- 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
  time0 = wrf_times_c(x, 0)
  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))
 

;---Read, vertical velocity(w), u and v winds, temp(tc), pressure(p), relative humidity(RH),one by one to reduce memory use
  
  p1 = wrf_user_getvar(a,"pressure",-1)   ; pressure
  p = dim_avg_n_Wrap(p1(startindx:endindx,:,:,:),0)
  delete(p1)

  tc1  = wrf_user_getvar(a,"tc",-1)    ; 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)

   
  rh1 = wrf_user_getvar(a,"rh",-1)        ; 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)
  
  u1   = wrf_user_getvar(a,"ua",-1)    ; 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)

  v1   = wrf_user_getvar(a,"va",-1)    ; 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)

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

;---------------------
      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       = "%"


     system("/bin/rm -f  file850wind.nc")
 
     ncdf = addfile("file850wind.nc" ,"c")  ; open output netCDF file


       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      = u_plane
       ncdf->v       = v_plane
       ncdf->w       = w_plane
     
end