load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

begin

; for use of high resolution map database, execute
; 'setenv NCARG_RANGS /network/rit/lab/atm419lab/zhou/rangs' (on snow)
; 'setenv NCARG_RANGS /atm419lab/zhou/rangs' (on ash, reed)
; before running the script

c="./cmap_RH.txt"
cmap=asciiread(c,(/72,3/),"float")
cmap=cmap/255

center=asciiread("../../../DATA_FAILURE_END/latlon_failure_end_sub.txt",(/53,2/),"float")
;print(center(:,1))
tinfo=asciiread("../../../DATA_FAILURE_END/time_failure_end_sub.txt",(/53,2/),"integer")
yyyymmddhh=tinfo(:,1)
;print(yyyymmddhh)
name=asciiread("../../../DATA_FAILURE_END/name_failure_sub.txt",(/53/),"string")
;print(name)

dcase=dimsizes(yyyymmddhh)
RH_comp=new((/dcase,161,161/),"float")

do case=0,dcase-1,1

print("working on "+yyyymmddhh(case)+" -- "+name(case))

data=addfile("../../../DATA_FAILURE_END/"+yyyymmddhh(case)+"_RH.grib","r")
;names=getfilevarnames(data)
;print(names)

RH:=data->R_GDS0_ISBL
RH!0="lev"
RH!1="lat"
RH!2="lon"
;printVarSummary(RH)
;print(num(ismissing(RH)))
;print(RH&lev)
;print(RH&lat)
;print(RH&lon)

lat=RH&lat
lon=RH&lon
centlat_ind=ind_nearest_coord(center(case,0),lat,0)
centlon_ind=ind_nearest_coord(center(case,1),lon,0)
centlat=lat(centlat_ind)
centlon=lon(centlon_ind)

dsizes=dimsizes(RH)
;print(dsizes)
dlat=dsizes(1)
dlon=dsizes(2)

RH_avg=dim_avg_n_Wrap(RH,0)

;printVarSummary(RH_avg)

RH_comp(case,:,:)=RH_avg(centlat_ind+80:centlat_ind-80,centlon_ind-80:centlon_ind+80)
;printVarSummary(RH_comp)

end do

rm=systemfunc("rm -rf RH_composite_sub.nc")

fout=addfile("./RH_composite_sub.nc","c")
fout->RH=RH_comp

end