;----------------------------------------------------------------------
; This plots all three data files using the WRF projection defined
; on the domain 01 file.
;----------------------------------------------------------------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

begin
;---Open all three files
files = "kdh_geo_em.d0" + (/1,2,3/) + ".nc"
a1 = addfile(files(0),"r")
a2 = addfile(files(1),"r")
a3 = addfile(files(2),"r")

;---Read variables off file into local arrays.

print("======First domain=====")
ter1 = a1->HGT_M(0,:,:)
lat1 = a1->XLAT_M(0,:,:)
lon1 = a1->XLONG_M(0,:,:)
printMinMax(ter1,0)
printMinMax(lat1,0)
printMinMax(lon1,0)


print("======Second domain=====")
ter2 = a2->HGT_M(0,:,:)
lat2 = a2->XLAT_M(0,:,:)
lon2 = a2->XLONG_M(0,:,:)
printMinMax(ter2,0)
printMinMax(lat2,0)
printMinMax(lon2,0)

print("======Third domain=====")
ter3 = a3->HGT_M(0,:,:)
lat3 = a3->XLAT_M(0,:,:)
lon3 = a3->XLONG_M(0,:,:)
printMinMax(ter3,0)
printMinMax(lat3,0)
printMinMax(lon3,0)

;----------------------------------------------------------------------
; Start the graphics
;----------------------------------------------------------------------
type = "png"
wks = gsn_open_wks(type,"panel_plot_geo_domains_gsn_native_dom1")

plot = new(3,graphic)

;---Set common resources for all three plots
res                      = True        ; Create some plot resources
res@gsnDraw              = False       ; Will draw plots later in panel
res@gsnFrame             = False
res@cnFillOn             = True        ; Create a color fill plot
res@cnFillPalette        = "rainbow"   ; Associate color map with contours

res@cnLinesOn            = False
res@cnLineLabelsOn       = False

res@cnLevelSelectionMode = "ManualLevels"
res@cnMinLevelValF       = 0.
res@cnMaxLevelValF       = 3500.
res@cnLevelSpacingF      = 250.

res@lbLabelBarOn         = False       ; Will draw later in panel
res@gsnAddCyclic         = False

res@mpOutlineBoundarySets = "AllBoundaries"  ; gives us more boundaries (for debuggig purposes)

;---Set the map projection using the projection info from the biggest domain (dom01)
res = wrf_map_resources(a1,res)
  
res@mpGeophysicalLineColor      = "Black"  ; Overwrite basic map settings
res@mpGridLineColor             = "Gray"
res@mpLimbLineColor             = "Gray"
res@mpNationalLineColor         = "Black"
res@mpPerimLineColor            = "Black"
res@mpUSStateLineColor          = "Black"
  
ter1@lat2d = lat1   ; This tells NCL where the data lies on the map
ter1@lon2d = lon1   ; This is necessary when you are not using
ter2@lat2d = lat2   ; a native projection.
ter2@lon2d = lon2
ter3@lat2d = lat3
ter3@lon2d = lon3

res@tiMainString  = files(0)
plot(0) = gsn_csm_contour_map(wks,ter1,res)

res@tiMainString  = files(1)
plot(1) = gsn_csm_contour_map(wks,ter2,res)

res@tiMainString  = files(2)
plot(2) = gsn_csm_contour_map(wks,ter3,res)

;---Set panel resources
pnlres                            = True
pnlres@gsnMaximize                = True
pnlres@gsnPanelYWhiteSpacePercent = 5       ; Add white space b/w plots.
pnlres@gsnPanelLabelBar           = True    ; Turn on common labelbar
pnlres@lbBoxMinorExtentF          = 0.15
pnlres@txString                   = "Using WRF projection from " + files(0)
pnlres@gsnPanelRowSpec            = True
gsn_panel(wks,plot,(/2,1/),pnlres)

end