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_START/latlon_failure_start_sub.txt",(/53,2/),"float")
;print(center(:,1))
tinfo=asciiread("../../../DATA_FAILURE_START/time_failure_start_sub.txt",(/53,2/),"integer")
yyyymmddhh=tinfo(:,1)
;print(yyyymmddhh)
name=asciiread("../../../DATA_FAILURE_START/name_failure_sub.txt",(/53/),"string")
;print(name)

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

do case=0,dcase-1,1

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

data=addfile("../../../DATA_FAILURE_START/ei.oper.an.pl.regn128sc."+yyyymmddhh(case),"r")
;names=getfilevarnames(data)
;print(names)

RH:=data->R_GDS4_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=new((/dlat,dlon/),"float")

do j=0,dlat-1,1
  do i=0,dlon-1,1
    RH_avg(j,i)=avg(RH(25:30,j,i)) ;{lev|25:30} corresponds to 700-850hPa
  end do
end do

copy_VarMeta(RH(0,:,:),RH_avg)
;printVarSummary(RH_avg)

RH_comp(case,:,:)=RH_avg(centlat_ind+21:centlat_ind-21,centlon_ind-21:centlon_ind+21)
;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