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/y98/julia/WRF-CLIM/Era-Int/wrfout_d01.nc","r") f2 = addfile("/scratch/y98/julia/WRF-CLIM/Era-Int/wrfout_d02.nc","r") f3 = addfile("/scratch/y98/julia/WRF-CLIM/Era-Int/wrfout_d03.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) ; get grid3 lat, lon wrf_lat2d3 = f3->XLAT(0,:,:) wrf_lon2d3 = f3->XLONG(0,:,:) topo3 = wrf_user_getvar(f3,"HGT",0) landmask3 = wrf_user_getvar(f3,"LANDMASK",0) topo3_mask = mask(topo3,landmask3,1) ;;convert lon2d to be 0->360 not -180->180 wrf_lon2d1 = where(wrf_lon2d1.lt.0,360.+wrf_lon2d1,wrf_lon2d1) wrf_lon2d2 = where(wrf_lon2d2.lt.0,360.+wrf_lon2d2,wrf_lon2d2) wrf_lon2d3 = where(wrf_lon2d3.lt.0,360.+wrf_lon2d3,wrf_lon2d3) dsizes_topo2 = dimsizes(topo2) nlat2 = dsizes_topo2(0) nlon2 = dsizes_topo2(1) dsizes_topo3 = dimsizes(topo3) nlat3 = dsizes_topo3(0) nlon3 = dsizes_topo3(1) ; get boundaries of grid 2 dsizes_grid2 = dimsizes(wrf_lat2d2) ;print(dsizes_grid2) s_lon2 = wrf_lon2d2(0,:) s_lat2 = wrf_lat2d2(0,:) n_lon2 = wrf_lon2d2(dsizes_grid2(0)-1,:) n_lat2 = wrf_lat2d2(dsizes_grid2(0)-1,:) w_lon2 = wrf_lon2d2(:,0) w_lat2 = wrf_lat2d2(:,0) e_lon2 = wrf_lon2d2(:,dsizes_grid2(1)-1) e_lat2 = wrf_lat2d2(:,dsizes_grid2(1)-1) ; get boundaries of grid 3 dsizes_grid3 = dimsizes(wrf_lat2d3) ;print(dsizes_grid3) s_lon3 = wrf_lon2d3(0,:) s_lat3 = wrf_lat2d3(0,:) n_lon3 = wrf_lon2d3(dsizes_grid3(0)-1,:) n_lat3 = wrf_lat2d3(dsizes_grid3(0)-1,:) w_lon3 = wrf_lon2d3(:,0) w_lat3 = wrf_lat2d3(:,0) e_lon3 = wrf_lon2d3(:,dsizes_grid3(1)-1) e_lat3 = wrf_lat2d3(:,dsizes_grid3(1)-1) dsizes_topo1 = dimsizes(topo1) nlat1 = dsizes_topo1(0) nlon1 = dsizes_topo1(1) pole_lat1 = 60.31 pole_lon1 = 321.38 type = "pdf" ;type = "x11" ;type@wkColorMap = "gui_default" ; define color scheme to use for contouring wks = gsn_open_wks(type,"Plot/grid_map") gsn_define_colormap(wks,"BlGrYeOrReVi200") res = True ; use color res@tiMainString = "(a)" 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 = "CylindricalEquidistant" res@mpCenterRotF = -3. res@mpCenterLatF = pole_lat1-90.0 res@mpCenterLonF = pole_lon1-180.0 res@mpOutlineDrawOrder = "PostDraw" ; force map tp be drawn 1st res@mpGridLineDashPattern = 2 ; lat/lon lines as dashed res@mpPerimOn = True res@mpPerimDrawOrder = "PostDraw" res@mpOutlineOn = True res@mpGridAndLimbOn = True res@mpGridSpacingF = 10. res@pmTickMarkDisplayMode = "Always" res@tmXTLabelsOn = False res@tmYRLabelsOn = False res@mpOutlineBoundarySets = "National" res@mpGeophysicalLineThicknessF = 1.5 res@mpDataBaseVersion = "MediumRes" res@mpDataSetName = "Earth..4" res@mpOutlineSpecifiers = "Australia:states" res@mpFillOn = True ;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 plres@gsLineLabelString= "test" dum = new(dimsizes(s_lon2),graphic) dum1 = new(dimsizes(n_lon2),graphic) dum2 = new(dimsizes(w_lon2),graphic) dum3 = new(dimsizes(e_lon2),graphic) dum = gsn_add_polyline(wks,map,s_lon2,s_lat2,plres) dum1 = gsn_add_polyline(wks,map,n_lon2,n_lat2,plres) dum2 = gsn_add_polyline(wks,map,w_lon2,w_lat2,plres) dum3 = gsn_add_polyline(wks,map,e_lon2,e_lat2,plres) polyres = True dum_3 = new(dimsizes(s_lon3),graphic) dum1_3 = new(dimsizes(n_lon3),graphic) dum2_3 = new(dimsizes(w_lon3),graphic) dum3_3 = new(dimsizes(e_lon3),graphic) dum_3 = gsn_add_polyline(wks,map,s_lon3,s_lat3,plres) dum1_3 = gsn_add_polyline(wks,map,n_lon3,n_lat3,plres) dum2_3 = gsn_add_polyline(wks,map,w_lon3,w_lat3,plres) dum3_3 = gsn_add_polyline(wks,map,e_lon3,e_lat3,plres) polyres = True polyres@gsMarkerIndex = 16 ;dum4 = gsn_add_polymarker(wks,map,115.8589,-31.9522,polyres) txres = True txres@txFontHeightF = 0.009 dum_t1 = gsn_add_text(wks,map,"Grid 2",113,-40.5,txres) dum_t2 = gsn_add_text(wks,map,"Grid 3",116.5,-36.5,txres) ;dum_t3 = gsn_add_text(wks,map,"(a)",97.5,-17.5,txres) ;draw(map) ;frame(wks) delete(res) 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_lat2d3(0,0) res@mpLeftCornerLonF = wrf_lon2d3(0,0) res@mpRightCornerLatF = wrf_lat2d3(nlat3-1,nlon3-1) res@mpRightCornerLonF = wrf_lon2d3(nlat3-1,nlon3-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 = "CylindricalEquidistant" res@mpCenterRotF = -3. res@mpCenterLatF = pole_lat1-90.0 res@mpCenterLonF = pole_lon1-180.0 res@mpOutlineDrawOrder = "PostDraw" ; force map tp be drawn 1st res@mpGridLineDashPattern = 2 ; lat/lon lines as dashed res@mpPerimOn = True res@mpPerimDrawOrder = "PostDraw" res@mpOutlineOn = True res@mpGridAndLimbOn = True res@mpGridSpacingF = 10. res@pmTickMarkDisplayMode = "Always" res@tmXTLabelsOn = False res@tmYRLabelsOn = False res@mpOutlineBoundarySets = "National" res@mpGeophysicalLineThicknessF = 1.5 res@mpDataBaseVersion = "MediumRes" res@mpDataSetName = "Earth..4" res@mpOutlineSpecifiers = "Australia:states" res@mpFillOn = True res@tiMainString = "(b)" 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 map2 = gsn_csm_contour_map(wks,topo3_mask,res) dum4 = gsn_add_polymarker(wks,map2,115.8589,-31.9522,polyres) dum_t3 = gsn_add_text(wks,map2,"Perth",116.2,-32,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/),(/2,1/),resP) ;gsn_panel(wks,(/map,map2,plot/),(/1,3/),resP) delete(wks) system("pdfcrop Plot/grid_map.pdf Plot/grid_map.pdf") end