load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "/home/users/zhangyl/fore_fnl/ncl_script/wrflib.ncl"
load "/home/users/zhangyl/fore_fnl/ncl_script/write_netcdf.ncl"
begin
; The WRF ARW input file.
  wrfdir = "/mnt/tianyuan/zhangyl/wrf_history/tuerqi/"
  extradir = "/mnt/tianyuan/zhangyl/extract_nc/tuerqi/10/"
  files = systemfunc("ls " + wrfdir + "wrfout*")
  nfile  = dimsizes(files)
;=================================================================================
; The specific height levels that we want the data interpolated to.
  height_levels = (/ 10 /)      ; height levels to plot - in meter
  height = height_levels
  nheight = dimsizes(height)
  height_levels@description = "Altitude"
  height_levels@units    = "m"
;==================================================================================
  filename = extradir+"2015_rh.nc"
  system("rm -f " + filename)
  fout = addfile(filename,"c")
;====================================================================================
 ;First get the variables we will need
  a1 = addfile(files(0),"r")
  lat   = wrf_user_getvar(a1,"lat",-1)
  lon   = wrf_user_getvar(a1,"lon",-1)
  dimz  = dimsizes(lat)
  nlat  = dimz(1)
  nlon  = dimz(2)
  ntime = 8760
  nlev  = 38
  z  = new((/ntime,nlev,nlat,nlon/),"float")
  rh = new((/ntime,nlev,nlat,nlon/),"float")
 do time=0,119,1
  rh(time,:,:,:)  = wrf_user_getvar(a1,"rh",time)        ; relative humidity
  z(time,:,:,:)  = wrf_user_getvar(a1,"z",time)
  z(time,:,:,:)  = z(time,:,:,:) - conform(z(time,:,:,:), wrf_user_getvar(a1, "ter", time), (/1, 2/)) 
 end do
;====================================================================================
 do i=1,nfile-1,1
  a = addfile(files(i),"r")
  do time=6,125,1
  rh(120*i+time-6,:,:,:) = wrf_user_getvar(a,"rh",time)
  z(120*i+time-6,:,:,:) = wrf_user_getvar(a,"z",time)
  z(120*i+time-6,:,:,:) = z(120*i+time-6,:,:,:) - conform(z(120*i+time-6,:,:,:), wrf_user_getvar(a, "ter", time), (/1, 2/))
  end do 
 end do
  rh_plane = new((/ntime,nheight,nlat,nlon/),"float")
  rh_plane  = wrf_user_intrp3d(rh,z,"h",height,0.,False)
  delete([/rh, z /])
  rh_plane@description = "Relative Humidity" ;
  rh_plane@units = "%"  ;
;====================================================================================
  fAtt               = True
  fAtt@title         = "NCL Simple Approach to netCDF Creation"
  fAtt@source_file   =  "a"
  fAtt@Conventions   = "None"
  fAtt@creation_date = systemfunc ("date")
  re1 = CopyGlobalAtts(fout,a1)
  re2 = CopyDimVar(fout,a1,height_levels)

 ;filevardef (fout,"Times",typeof(Times),getvardims(Times))
  filevardef (fout,"Height",typeof(height_levels),(/"height"/))
  filevardef (fout, "XLAT" ,typeof(lat),getvardims(lat))
  filevardef (fout, "XLONG" ,typeof(lon),getvardims(lon))

  varvar_types2D = (/ "float", "float",  "float" /)
  var_names2D = (/ "Height", "XLAT", "XLONG" /)
  
  varvar_types4D = (/ "float" /)
  var_names4D = (/ "RH" /)
  filevardef (fout, var_names4D , varvar_types4D , (/"Time","height","south_north","west_east"/) )
 
 ;filevarattdef(fout,"Times",Times)                    ; copy Times attributes
  filevarattdef(fout,"Height",height_levels)           ; copy height attributes
  filevarattdef(fout,"XLAT",lat)                       ; copy lat attributes
  filevarattdef(fout,"XLONG",lon)                      ; copy lon attributes
  filevarattdef(fout,"RH",rh_plane)                    ; copy rh_plane attributes
;===============================================================================
 ;fout->Times   = Times
  fout->Height  = height_levels
  fout->XLAT    = (/lat/)
  fout->XLONG   = (/lon/)
  fout->RH    = (/rh_plane/)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  delete(rh_plane)
end