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/wrf/WRFUserARW.ncl"

begin
;f1 = addfile("/scratch/partner69/jatinkala/WRF3.3/WRFV3/run/test-fnl-rrtmg-mp6/wrfout_d01_2009-10-07_00:00:00.nc","r")
;f2 = addfile("/scratch/partner69/jatinkala/WRF3.3/WRFV3/run/test-fnl-rrtmg-mp6/wrfout_d02_2009-10-07_00:00:00.nc","r")
f1 = addfile("/scratch/y98/jatinkala/WRF3.3/WRFV3/run/test-nnrp/wrfinput_d01" + ".nc","r")
f2 = addfile("/scratch/y98/jatinkala/WRF3.3/WRFV3/run/test-nnrp/wrfout_d02_2010-01-01_00:00:00" + ".nc","r")

; get grid 1 topo, land maks, lat, lon
topo1 = wrf_user_getvar(f1,"HGT",0)
landmask1 = wrf_user_getvar(f1,"LANDMASK",0)
wrf_lat2d1 = f1->XLAT(0,:,:)
wrf_lon2d1 = f1->XLONG(0,:,:)
topo1_mask = mask(topo1,landmask1,1)

; get grid2 lat, lon
wrf_lat2d2 = f2->XLAT(0,:,:)
wrf_lon2d2 = f2->XLONG(0,:,:)
topo2 = wrf_user_getvar(f2,"HGT",0)
landmask2 = wrf_user_getvar(f2,"LANDMASK",0)
topo2_mask = mask(topo2,landmask2,1)

dsizes_topo2 = dimsizes(topo2)
nlat2 = dsizes_topo2(0)
nlon2 = dsizes_topo2(1)

; get boundaries of grid 2
dsizes_grid2 = dimsizes(wrf_lat2d2)
;print(dsizes_grid2)
s_lon = wrf_lon2d2(0,:)
s_lat = wrf_lat2d2(0,:)
n_lon = wrf_lon2d2(dsizes_grid2(0)-1,:)
n_lat = wrf_lat2d2(dsizes_grid2(0)-1,:)
w_lon = wrf_lon2d2(:,0)
w_lat = wrf_lat2d2(:,0)
e_lon = wrf_lon2d2(:,dsizes_grid2(1)-1)
e_lat = wrf_lat2d2(:,dsizes_grid2(1)-1)

dsizes_topo1 = dimsizes(topo1)
nlat1 = dsizes_topo1(0)
nlon1 = dsizes_topo1(1)

type = "pdf"
;type = "x11"
;type@wkColorMap = "gui_default" ; define color scheme to use for contouring
wks = gsn_open_wks(type,"eps-figs/grid_map")

res = True ; use color
res@cnLinesOn = False ; no lines on contour plot
res@cnFillOn = True ; fill map with color
res@cnLineLabelsOn = False ; no line label
res@gsnSpreadColors = True ; use full spread of color-bar
; define corners of plot, since we are not plotting whole grid (not plotting over ocean)
res@mpLimitMode       = "Corners"
res@mpLeftCornerLatF  = wrf_lat2d1(0,0)
res@mpLeftCornerLonF  = wrf_lon2d1(0,0)
res@mpRightCornerLatF = wrf_lat2d1(nlat1-1,nlon1-1)
res@mpRightCornerLonF = wrf_lon2d1(nlat1-1,nlon1-1)
; set some more basic plot aspects
res@mpGridLatSpacingF = 10 ; grid line spacing on plot
res@mpGridLonSpacingF = 10
res@pmTickMarkDisplayMode = "Always" ; display tick marks
res@tmXTLabelsOn = False ; turn off tick-marks labels at top of plot
res@tmYRLabelsOn = False ; trun off tick-marks labels at right of plot
res@tmXTOn = False ; turn off top tick marks
res@tmYROn = False; turn off right tick marks
res@mpGridAndLimbOn = True ; set to true if you want grid lines to show
res@mpGridLineDashPattern = 2 ; make grid lines dash
res@mpGridLineColor = "black" ; grid line color black
res@mpGeophysicalLineThicknessF = 1.5
res@cnLevelSelectionMode = "ManualLevels"
res@cnMaxLevelValF = 800.
res@cnMinLevelValF = 50.
res@cnLevelSpacingF = 50.
res@lbLabelBarOn = False

; define map projection params
truelat1 = f1@TRUELAT1
truelat2 = f1@TRUELAT2
stand_lon = f1@CEN_LON
res@mpProjection        = "LambertConformal"
res@mpLambertParallel1F = truelat1
res@mpLambertParallel2F = truelat2
res@mpLambertMeridianF  = stand_lon
res@tfDoNDCOverlay = True ; this is necessary since grid is already native

res@gsnDraw = False
res@gsnFrame = False

res@lbTitleString = "m"
res@lbTitlePosition = "Bottom"

;type@wkColorMap = "gui_default" ; define color scheme to use for contouring

map = gsn_csm_contour_map(wks,topo1_mask,res)
plres = True
plres@gsLineThicknessF = 2.5
dum = new(dimsizes(s_lon),graphic)
dum1 = new(dimsizes(n_lon),graphic)
dum2 = new(dimsizes(w_lon),graphic)
dum3 = new(dimsizes(e_lon),graphic)

dum = gsn_add_polyline(wks,map,s_lon,s_lat,plres)
dum1 = gsn_add_polyline(wks,map,n_lon,n_lat,plres)
dum2 = gsn_add_polyline(wks,map,w_lon,w_lat,plres)
dum3 = gsn_add_polyline(wks,map,e_lon,e_lat,plres)
polyres = True
polyres@gsMarkerIndex = 16
dum4 = gsn_add_polymarker(wks,map,115.8589,-31.9522,polyres)

txres = True
txres@txFontHeightF = 0.023
dum_t1 = gsn_add_text(wks,map,"Grid 2",117.5,-41,txres)
dum_t2 = gsn_add_text(wks,map,"Perth",111.8,-31.9522,txres)
dum_t3 = gsn_add_text(wks,map,"(a)",97.5,-17.5,txres)

;draw(map)
;frame(wks)

delete(res)

; plot topo2
res = True
res@cnLinesOn = False ; no lines on contour plot
res@cnFillOn = True ; fill map with color
res@cnLineLabelsOn = False ; no line label
res@gsnSpreadColors = True ; use full spread of color-bar
; define corners of plot, since we are not plotting whole grid (not plotting over ocean)
res@mpLimitMode       = "Corners"
res@mpLeftCornerLatF  = wrf_lat2d2(0,0)
res@mpLeftCornerLonF  = wrf_lon2d2(0,0)
res@mpRightCornerLatF = wrf_lat2d2(nlat2-1,nlon2-1)
res@mpRightCornerLonF = wrf_lon2d2(nlat2-1,nlon2-1)
; set some more basic plot aspects
res@mpGridLatSpacingF = 5 ; grid line spacing on plot
res@mpGridLonSpacingF = 5
res@pmTickMarkDisplayMode = "Always" ; display tick marks
res@tmXTLabelsOn = False ; turn off tick-marks labels at top of plot
res@tmYRLabelsOn = False ; trun off tick-marks labels at right of plot
res@tmXTOn = False ; turn off top tick marks
res@tmYROn = False; turn off right tick marks
res@mpGridAndLimbOn = True ; set to true if you want grid lines to show
res@mpGridLineDashPattern = 2 ; make grid lines dash
res@mpGridLineColor = "black" ; grid line color black
res@mpGeophysicalLineThicknessF = 1.5

; define map projection params
truelat12 = f2@TRUELAT1
truelat22= f2@TRUELAT2
stand_lon2 = f2@CEN_LON
res@mpProjection        = "LambertConformal"
res@mpLambertParallel1F = truelat12
res@mpLambertParallel2F = truelat22
res@mpLambertMeridianF  = stand_lon2
res@tfDoNDCOverlay = True ; this is necessary since grid is already native

res@gsnDraw = False
res@gsnFrame = False

res@lbTitleString = "m"
res@lbTitlePosition = "Bottom"

res@cnLevelSelectionMode = "ManualLevels"
res@cnMaxLevelValF = 800.
res@cnMinLevelValF = 50.
res@cnLevelSpacingF = 50.
res@lbLabelBarOn = False

;type@wkColorMap = "gui_default" ; define color scheme to use for contouring

map2 = gsn_csm_contour_map(wks,topo2_mask,res)
dum_t4 = gsn_add_text(wks,map2,"(b)",109.,-28.,txres)
dum_DS = gsn_add_text(wks,map2,"Darling~C~scarp",112.5,-32,txres)

delete(res)

;Darling Scarp plot
f = addfile("/group/y98/jatinkala/GEO-SCI-AUS-TOPO/aus_dem_9s.nc","r")
topo = f->height
lat = f->latitude
lon = f->longitude
ind_lat = ind(lat .ge. -34. .and. lat .le. -31.)
ind_lon = ind(lon .ge. 114.5 .and. lon .le. 117.0)
scarp_topo = topo(ind_lat,ind_lon)
scarp_real_topo = scarp_topo*scarp_topo@scale_factor + scarp_topo@add_offset
delete(scarp_topo)
scarp_lat = lat(ind_lat)
scarp_lon = lon(ind_lon)
lat2d = conform_dims(dimsizes(scarp_real_topo),scarp_lat,0)
lon2d = conform_dims(dimsizes(scarp_real_topo),scarp_lon,1)
scarp_real_topo@lat2d = lat2d
scarp_real_topo@lon2d = lon2d

; add boundaries of darling scarp onto map2 here
xdum = new(2,graphic)
xdum1 = new(2,graphic)
xdum2 = new(2,graphic)
xdum3 = new(2,graphic)
xdum = gsn_add_polyline (wks,map2,(/min(scarp_lon),max(scarp_lon)/),(/min(scarp_lat),min(scarp_lat)/),plres)
xdum1 = gsn_add_polyline(wks,map2,(/min(scarp_lon),max(scarp_lon)/),(/max(scarp_lat),max(scarp_lat)/),plres)
xdum2 = gsn_add_polyline(wks,map2,(/min(scarp_lon),min(scarp_lon)/),(/min(scarp_lat),max(scarp_lat)/),plres)
xdum3 = gsn_add_polyline(wks,map2,(/max(scarp_lon),max(scarp_lon)/),(/min(scarp_lat),max(scarp_lat)/),plres)

; plot here
res = True
res@gsnDraw = False
res@gsnFrame = False
res@cnLinesOn = False
res@cnFillOn = True
res@cnLineLabelsOn = False
res@mpGridAndLimbOn        = True
res@pmTickMarkDisplayMode  = "Always"           ; turn on tickmarks
res@mpGridSpacingF = 0.5
res@mpGridLineDashPattern = 2 ; make grid lines dash
res@mpGridLineColor = "black" ; grid line color black
res@mpMaxLatF = max(lat2d)
res@mpMinLatF = min(lat2d)
res@mpMaxLonF = max(lon2d)
res@mpMinLonF = min(lon2d)
res@mpDataBaseVersion    = "HighRes"
res@cnLevelSelectionMode = "ManualLevels"
res@cnMaxLevelValF = 750.
res@cnMinLevelValF = 50.
res@cnLevelSpacingF = 50.
res@lbLabelBarOn = False

;wks = gsn_open_wks("eps","eps-figs/Darling_Scarp")
plot = gsn_csm_contour_map_ce(wks,scarp_real_topo,res)
txres@txFontHeightF = 0.02
dum_t5 = gsn_add_text(wks,plot,"(C)",114.75,-31.25,txres)

resP = True
resP@gsnPanelLabelBar = True
resP@lbTitleString = "(m)"
resP@lbTitlePosition = "Bottom"
resP@lbTitleFontHeightF = 0.02
resP@pmLabelBarOrthogonalPosF = -0.015

gsn_panel(wks,(/map,map2,plot/),(/1,3/),resP)
delete(wks)
system("pdfcrop eps-figs/grid_map.pdf eps-figs/grid_map.pdf")

end