;; ==============================================
;; bug_wrf_user_ll_to_ij_and_rcm2points.ncl
;;
;; Written by: Jared A. Lee
;;             jaredlee@ucar.edu
;; Written on: 21 Nov 2017
;;
;; This script demonstrates a likely bug in the wrf_user_ll_to_ij and rcm2points functions.
;; wrf_user_ll_to_ij returns non-existent (i,j) points when fed a (lon,lat) outside the WRF domain.
;; rcm2points(_Wrap) interpolates WRF model data to points WAY outside the model domain.
;; ==============================================

begin

wrf_dir		= "/glade/p/nral0011/ALASKA/WRFV3.7.1/PRODUCE_4KM_V2"
obs_dir		= "/glade/p/ral/RHAP/anewman/madis/2002/09/01"
wrf_fname	= wrf_dir+"/surface_d01_2002-09-01_00:00:00.nc"
obs_fname	= obs_dir+"/20020901_0000_metar"

;; Read in WRF data
print("Reading WRF file "+wrf_fname)
wrf_file	= addfile(wrf_fname, "r")

wrf_lat	= wrf_user_getvar(wrf_file, "XLAT", 0)
wrf_lon	= wrf_user_getvar(wrf_file, "XLONG", 0)
wrf_uv10	= wrf_user_getvar(wrf_file, "uvmet10", -1)
wrf_u10	= wrf_uv10(0,0,:,:)
wrf_v10	= wrf_uv10(1,0,:,:)
wrf_ws10	= sqrt(wrf_u10^2 + wrf_v10^2)
copy_VarMeta(wrf_u10, wrf_ws10)
wrf_ws10@description	= "10-m wind speed rotated to earth coordinates"
delete([/wrf_uv10, wrf_u10, wrf_v10/])

wrf_lat_min	= min(wrf_lat)
wrf_lat_max	= max(wrf_lat)
wrf_lon_min	= min(wrf_lon)
wrf_lon_max	= max(wrf_lon)
print("   wrf_lat_min = "+wrf_lat_min+" N, wrf_lat_max = "+wrf_lat_max+" N")
print("   wrf_lon_min = "+wrf_lon_min+" E, wrf_lon_max = "+wrf_lon_max+" E")

wrf_dims	= getfiledimsizes(wrf_file)
wrf_n_lon	= wrf_dims(2)
wrf_n_lat	= wrf_dims(3)
print("   wrf_n_lon = "+wrf_n_lon+", wrf_n_lat = "+wrf_n_lat)


;; Read in MADIS obs data
zip_obs	= False
print("Reading obs file "+obs_fname)
if (fileexists(obs_fname+".gz")) then
	zip_obs	= True
	system("gunzip "+obs_fname+".gz")
end if
obs_file	= addfile(obs_fname, "r")

obs_all_stn_name	= chartostring(obs_file->stationName)
obs_all_loc_name	= chartostring(obs_file->locationName)
obs_all_lat			=              obs_file->latitude
obs_all_lon			=              obs_file->longitude
obs_all_wspd		=              obs_file->windSpeed

;; wrf_user_ll_to_ij gives nonsense (i,j) values for points outside the domain
obs_all_ij	= wrf_user_ll_to_ij(wrf_file, obs_all_lon, obs_all_lat, True)

;; rcm2points interpolates to many points WAY outside the domain
wrf_ws10_obs_pts	= rcm2points_Wrap(wrf_lat, wrf_lon, wrf_ws10, obs_all_lat, obs_all_lon, 0)
wrf_ws10_obs_pts@_FillValue = -999

;; Show some of the output of these functions
print("")
print("All obs locations:")
print(obs_all_stn_name+", (Lon,Lat) = ("+obs_all_lon+" E, "+obs_all_lat+" N), Nearest (i,j) = ("+obs_all_ij(0,:)+", "+obs_all_ij(1,:)+")"+\
	", Interpolated WS10 = "+wrf_ws10_obs_pts+" m/s, Station Name = "+obs_all_loc_name)

;; Screen out the points where the interpolated value is missing
inds1	= ind(.not.ismissing(wrf_ws10_obs_pts))

wrf_ws10_obs_pts_sub1	= wrf_ws10_obs_pts(inds1)
obs_sub1_stn_name	= obs_all_stn_name(inds1)
obs_sub1_loc_name	= obs_all_loc_name(inds1)
obs_sub1_lat		= obs_all_lat(inds1)
obs_sub1_lon		= obs_all_lon(inds1)
obs_sub1_wspd		= obs_all_wspd(inds1)
obs_sub1_ij	= wrf_user_ll_to_ij(wrf_file, obs_sub1_lon, obs_sub1_lat, True)

print("")
print("Obs locations without missing interpolated WS10 from rcm2points:")
print(obs_sub1_stn_name+", (Lon,Lat) = ("+obs_sub1_lon+" E, "+obs_sub1_lat+" N), Nearest (i,j) = ("+obs_sub1_ij(0,:)+", "+obs_sub1_ij(1,:)+")"+\
	", Interpolated WS10 = "+wrf_ws10_obs_pts_sub1+" m/s, Station Name = "+obs_sub1_loc_name)

;; Screen out points with (i,j) within the bounds of the domain
inds2	= ind(obs_sub1_ij(0,:).lt.0 .or. obs_sub1_ij(0,:).gt.wrf_n_lon .or. obs_sub1_ij(1,:).lt.0 .or. obs_sub1_ij(1,:).gt.wrf_n_lat)

wrf_ws10_obs_pts_sub2	= wrf_ws10_obs_pts_sub1(inds2)
obs_sub2_stn_name	= obs_sub1_stn_name(inds2)
obs_sub2_loc_name	= obs_sub1_loc_name(inds2)
obs_sub2_lat		= obs_sub1_lat(inds2)
obs_sub2_lon		= obs_sub1_lon(inds2)
obs_sub2_wspd		= obs_sub1_wspd(inds2)
obs_sub2_ij	= wrf_user_ll_to_ij(wrf_file, obs_sub2_lon, obs_sub2_lat, True)

print("")
print("Obs locations with nearest (i,j) values outside the bounds of the WRF domain:")
print(obs_sub2_stn_name+", (Lon,Lat) = ("+obs_sub2_lon+" E, "+obs_sub2_lat+" N), Nearest (i,j) = ("+obs_sub2_ij(0,:)+", "+obs_sub2_ij(1,:)+")"+\
	", Interpolated WS10 = "+wrf_ws10_obs_pts_sub2+" m/s, Station Name = "+obs_sub2_loc_name)

print("")
print("Recall that wrf_n_lon = "+wrf_n_lon+", wrf_n_lat = "+wrf_n_lat)

if (zip_obs.eq.True) then
	system("gzip "+obs_fname)
end if

end