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" load "/scratch/y98/shared/ncl-functions/cordex_grid_plot_functions.ncl" begin x = 100 y = 100 domain = "d03" files_bom = systemfunc("ls /group/y98/julia/BOM/Processed/BOM_clim_means_" + domain + ".nc") f_bom = addfile(files_bom,"r") tmax_bom = f_bom->TMAX_mean_season tmin_bom = f_bom->TMIN_mean_season ppt_bom = f_bom->ppt_mean_season delete(f_bom) ;get a landmask and some other stuff from a wrfout file ; load a wrfout file to get land-mask info out f = addfile("/scratch/y98/julia/WRF-CLIM/Era-Int/wrfout_" + domain + ".nc","r") landmask = wrf_user_getvar(f,"LANDMASK",0) ;dsizes_array = dimsizes(land_mask) ; also get some map info for later on truelat1 = f@TRUELAT1 truelat2 = f@TRUELAT2 cen_lon = f@CEN_LON cen_lat = f@CEN_LAT wrf_lat2d = f->XLAT(0,:,:) wrf_lon2d = f->XLONG(0,:,:) latr = wrf_lat2d lonr = wrf_lon2d ;mask water and remove relaxation zone dsizes_array = dimsizes(landmask) print(dsizes_array) landmask_nr = landmask(10:(dsizes_array(0)-11),10:(dsizes_array(1)-11)) lon = lonr(10:(dsizes_array(0)-11),10:(dsizes_array(1)-11)) lat = latr(10:(dsizes_array(0)-11),10:(dsizes_array(1)-11)) dsizes_landmask_nr = dimsizes(landmask_nr) print(dsizes_landmask_nr) nlat = dsizes_landmask_nr(0) nlon = dsizes_landmask_nr(1) ; convert lon to be 0->360 not -180->180 lon = where(lon.lt.0,360.+lon,lon) pole_lat1 = 60.31 pole_lon1 = 321.38 tmax_bom_mask = mask(tmax_bom,landmask_nr,1) tmin_bom_mask = mask(tmin_bom,landmask_nr,1) ppt_bom_mask = mask(ppt_bom,landmask_nr,1) ;********* load bom station locations *********** f_bom = addfile("/scratch/y98/jatinkala/WRF3.3/WRFV3/run/bom-station-locs/bom_ppt_temp_ll.nc","r") ppt_lats1 = f_bom->ppt_lats1 ppt_longs1 = f_bom->ppt_longs1 tmax_lats1 = f_bom->tmax_lats1 tmax_longs1 = f_bom->tmax_longs1 ;*******Draw Panel Plots************ type= "pdf" wks_tmax = gsn_open_wks(type,"Plot/bom_tmax_" + domain) wks_tmin = gsn_open_wks(type,"Plot/bom_tmin_" + domain) wks_ppt = gsn_open_wks(type,"Plot/bom_ppt_" + domain) gsn_define_colormap(wks_tmax,"BlGrYeOrReVi200") gsn_define_colormap(wks_tmin,"BlGrYeOrReVi200") gsn_define_colormap(wks_ppt,"BlGrYeOrReVi200") ;TMAX** res = True ;define projection parameters 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 ;some more basic plot aspects 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 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@tmXBOn = False res@tmXTOn = False res@tmYLOn = False res@tmYROn = False 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 = 34. res@cnMinLevelValF = 12. res@cnLevelSpacingF = 1. res@lbLabelBarOn = False res@tfDoNDCOverlay = True ; this is necessary since grid is already native ; define corners of plot, since we are not plotting whole grid (not plotting over ocean) res@mpLimitMode = "Corners" res@mpLeftCornerLatF = lat(0,0) res@mpLeftCornerLonF = lon(0,0) res@mpRightCornerLatF = lat(nlat-1,nlon-1) res@mpRightCornerLonF = lon(nlat-1,nlon-1) res@gsnDraw = False res@gsnFrame = False res@lbTitleString = "C" res@lbTitlePosition = "Bottom" ;plot TMAX plot_tmax = new(4,graphic) plot_tmax(0) = gsn_csm_contour_map(wks_tmax,tmax_bom_mask(0,:,:),res) plot_tmax(1) = gsn_csm_contour_map(wks_tmax,tmax_bom_mask(1,:,:),res) plot_tmax(2) = gsn_csm_contour_map(wks_tmax,tmax_bom_mask(2,:,:),res) plot_tmax(3) = gsn_csm_contour_map(wks_tmax,tmax_bom_mask(3,:,:),res) ;add station locations ; add station locations dum = new(3000,graphic) counter = 0 pmres = True pmres@gsMarkerColor = "white" pmres@gsMarkerIndex = "dot" pmres@gsMarkerSizeF = 0.004 ; add ppt station locations, first check if stations within grid do i = 0,dimsizes(tmax_lats1)-1,1 if (tmax_lats1(i) .le. max(lat)) .and. (tmax_lats1(i) .ge. min(lat)) .and. (tmax_longs1(i) .le. max(lon)) .and. (tmax_longs1(i) .ge. min(lon)) dum(counter) = gsn_add_polymarker(wks_tmax,plot_tmax(0),tmax_longs1(i),tmax_lats1(i),pmres) counter = counter + 1 else end if end do ; now put 4 plots in panel resP = True resP@gsnPanelLabelBar = True ; put panel lable bar on (don't use this if still trying to figure out scale) resP@txString = "(a) Maximum Temperature" ; panel figure title ;resP@txFontHeightF = 0.035 resP@gsnPanelFigureStrings = (/"DJF","MAM","JJA","SON"/) ; text-box to go on plots showing season resP@gsnPanelFigureStringsFontHeightF = 0.014 resP@lbTitleString = "(~S~o~N~C)" resP@lbTitlePosition = "Bottom" resP@pmLabelBarOrthogonalPosF = -0.025 resP@lbTitleFontHeightF = 0.02 ;resP@pmLabelBarWidthF = 1.0 ;resP@amJust = "BottomLeft" ; put text-box bottom left gsn_panel(wks_tmax,plot_tmax,(/1,4/),resP) ; make panel plot system("pdfcrop Plot/bom_tmax_d03.pdf Plot/bom_tmax_d03.pdf") delete(resP) delete(dum) ;**TMIN** ; define label bar properties res@lbLabelBarOn = False ; Turn off label bar for individual plot, NOTE: set this to true if you are trying to figure out the scale ; set colorbar manually res@cnLevelSelectionMode = "ManualLevels" res@cnMinLevelValF = 2 res@cnMaxLevelValF = 16 res@cnLevelSpacingF = 1 plot_tmin = new(4,graphic) plot_tmin(0) = gsn_csm_contour_map(wks_tmin,tmin_bom_mask(0,:,:),res) plot_tmin(1) = gsn_csm_contour_map(wks_tmin,tmin_bom_mask(1,:,:),res) plot_tmin(2) = gsn_csm_contour_map(wks_tmin,tmin_bom_mask(2,:,:),res) plot_tmin(3) = gsn_csm_contour_map(wks_tmin,tmin_bom_mask(3,:,:),res) ; now put 4 plots in panel resP = True resP@gsnPanelLabelBar = True ; put panel lable bar on (don't use this if still trying to figure out scale) resP@txString = "(b) Minimum Temperature" ; panel figure title ;resP@txFontHeightF = 0.035 resP@gsnPanelFigureStrings = (/"DJF","MAM","JJA","SON"/) ; text-box to go on plots showing season resP@gsnPanelFigureStringsFontHeightF = 0.014 resP@lbTitleString = "(~S~o~N~C)" resP@lbTitlePosition = "Bottom" resP@pmLabelBarOrthogonalPosF = -0.025 resP@lbTitleFontHeightF = 0.02 ;resP@pmLabelBarWidthF = 1.0 ;resP@amJust = "BottomLeft" ; put text-box bottom left gsn_panel(wks_tmin,plot_tmin,(/1,4/),resP) ; make panel plot system("pdfcrop Plot/bom_tmin_d03.pdf Plot/bom_tmin_d03.pdf") delete(resP) ;**PPT** ; define label bar properties res@lbLabelBarOn = False ; Turn off label bar for individual plot, NOTE: set this to true if you are trying to figure out the scale ; set colorbar manually res@cnLevelSelectionMode = "ManualLevels" res@cnMinLevelValF = 20 res@cnMaxLevelValF = 200 res@cnLevelSpacingF = 20 plot_ppt = new(4,graphic) plot_ppt(0) = gsn_csm_contour_map(wks_ppt,ppt_bom_mask(0,:,:),res) plot_ppt(1) = gsn_csm_contour_map(wks_ppt,ppt_bom_mask(1,:,:),res) plot_ppt(2) = gsn_csm_contour_map(wks_ppt,ppt_bom_mask(2,:,:),res) plot_ppt(3) = gsn_csm_contour_map(wks_ppt,ppt_bom_mask(3,:,:),res) ; add station locations dum = new(3000,graphic) counter = 0 pmres = True pmres@gsMarkerColor = "white" pmres@gsMarkerIndex = "dot" pmres@gsMarkerSizeF = 0.004 ; add ppt station locations, first check if stations within grid do i = 0,dimsizes(ppt_lats1)-1,1 if (ppt_lats1(i) .le. max(lat)) .and. (ppt_lats1(i) .ge. min(lat)) .and. (ppt_longs1(i) .le. max(lon)) .and. (ppt_longs1(i) .ge. min(lon)) dum(counter) = gsn_add_polymarker(wks_ppt,plot_ppt(0),ppt_longs1(i),ppt_lats1(i),pmres) counter = counter + 1 else end if end do ; now put 4 plots in panel resP = True resP@gsnPanelLabelBar = True ; put panel lable bar on (don't use this if still trying to figure out scale) resP@txString = "(c) Rainfall" ; panel figure title ;resP@txFontHeightF = 0.035 resP@gsnPanelFigureStrings = (/"DJF","MAM","JJA","SON"/) ; text-box to go on plots showing season resP@gsnPanelFigureStringsFontHeightF = 0.014 resP@lbTitleString = "(mm month~S~-1~N~)" resP@lbTitlePosition = "Bottom" resP@pmLabelBarOrthogonalPosF = -0.025 resP@lbTitleFontHeightF = 0.02 ;resP@pmLabelBarWidthF = 1.0 ;resP@amJust = "BottomLeft" ; put text-box bottom left gsn_panel(wks_ppt,plot_ppt,(/1,4/),resP) ; make panel plot system("pdfcrop Plot/bom_ppt_d03.pdf Plot/bom_ppt_d03_crop.pdf") end