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/csm/contributed.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl" ;====================================================================== ; The main code ;====================================================================== begin ;---------------------------------------------------------- ; List Prism, WRF Present-Day, and WRF Deforested Directories ;---------------------------------------------------------- ;---Specify interpolation to be used method = "bilinear" wgtFileDir = "/glade/p/work/burakows/plot/NCL/ESMF_regridding/" wgtFileName = "PRISM_4kmM2_to_WRF."+method+"_wgt.nc" wgtFilePath = wgtFileDir+wgtFileName print("wgtFilePath="+wgtFilePath) ;---PRISM files [input] varp = "tmin" ; PRISM name dirp = "/glade/p/work/burakows/prism/data/mPRISM_"+varp+"_201111_201204/" filp = systemfunc("cd "+dirp+" ; ls PRISM_"+varp+"*bil.nc") nfilp = dimsizes(filp) print(filp(0:nfilp-1)) ; 1st 5 file ;---WRF Deforested directory varw = "T2min" ; WRF variable name dirDF = "/glade/p/work/burakows/plot/LC-1850/" dirDFs = systemfunc("cd "+dirDF+" ; ls -d albedo_MP*_201111-201204") ;---WRF Present-Day directory dirPD = "/glade/p/work/burakows/plot/Present-Day/" dirPDs = systemfunc("cd "+dirPD+" ; ls -d albedo_MP*_201111-201204") ndirPD = dimsizes(dirPDs) ; number of wrf directories prtFlg = True wrfPD = addfile(dirPD+dirPDs(0)+"/"+"wrfout_d03_2011-11-01_00:00:00.nc","r") ; read any wrfout file x = wrfPD->T2 ; read T2 and preserve metadata x@lat2d = wrfPD->XLAT(0,:,:) ; Assign lat to x x@lon2d = wrfPD->XLONG(0,:,:) ; Assign lon to x minlat = min(x@lat2d) maxlat = max(x@lat2d) minlon = min(x@lon2d) maxlon = max(x@lon2d) bndadd = 0.25 ; spacing around the plot edges ;--- Read in wrf land mask l = wrfPD->LANDMASK lm = l(0,:,:) wrfDF = addfile(dirDF+dirDFs(0)+"/"+"wrfout_d03_2011-11-01_00:00:00.nc","r") ; read any wrfout file ;---------------------------------------------------------------------- ; Create Land Cover Type Conversion masks ;---------------------------------------------------------------------- luPD = wrfPD->LU_INDEX(0,:,:) ; Present-Day Land Cover luPD@_FillValue = -9999 ; Assign fill value luDF = wrfDF->LU_INDEX(0,:,:) ; 1850 Deforested Land Cover ludiff = luPD ; copy metadata to difference in land cover ludiff = luPD-luDF ; Calculate difference in land cover ludiff@lat2d = wrfPD->XLAT(0,:,:) ludiff@lon2d = wrfPD->XLONG(0,:,:) ludiff = where(ludiff.eq.10 .or. \ ludiff.eq.9 .or. \ ludiff.eq.6 .or. \ ludiff.eq.-4.or. \ ludiff.eq.1., ludiff, ludiff@_FillValue) printVarSummary(ludiff) ;---------------------------------------------------------------------- ; Plotting options section ;---------------------------------------------------------------------- pltType = "ps" ; plot type pltDir = "./" ; plot directory pltName = "Panel_biascorr_"+varw+"PD-DF_2011" ; plot name (ps file) pltPath = pltDir+pltName ; plot path wks = gsn_open_wks(pltType,pltPath) ; create workstation for ps file gsn_define_colormap(wks,"temp_diff_18lev") res = True res@gsnDraw = False res@gsnFrame = False res@cnFillOn = True ; color plot desired res@cnLinesOn = False ; turn off contour lines res@cnLineLabelsOn = False ; turn off contour labels res@cnInfoLabelOn = False ; turn off info label (top labels of indvid. plots) res@cnFillMode = "RasterFill" ; turn raster on res@cnLevelSelectionMode = "ManualLevels" ; Set contour levels manually res@cnMinLevelValF = -4 ; minimum contour (degrees C, or mm) res@cnMaxLevelValF = 4 ; maximum contour (degrees C, or mm) res@cnLevelSpacingF = 0.5 ; contour interval res@lbLabelBarOn = False ; Will turn on in panel later res@lbOrientation = "Horizontal" ; Horizontal label bar res@mpFillOn = False res@mpOutlineOn = True res@mpOutlineBoundarySets = "AllBoundaries" res@mpProjection = "CylindricalEquidistant" res@mpLimitMode = "LatLon" ; required res@mpMinLatF = minlat-bndadd res@mpMaxLatF = maxlat+bndadd res@mpMinLonF = minlon-bndadd res@mpMaxLonF = maxlon+bndadd res@mpCenterLonF = (minlon + maxlon)*0.5 res@mpCenterLatF = (minlat + maxlat)*0.5 res@gsnLeftString = "" ; Turn off left subtitle res@gsnRightString = "" ; Turn off right subtitle res@gsnMajorLatSpacing = 2 ;res@gsnMajorLonSpacing = 2 res@gsnMinorLonSpacing = 2 res@gsnAddCyclic = False ; regional grid (changes central meridian)/xwo ;----Resources for stipling of significant signal/noise (>2) snres = True snres@gsnDraw = False snres@gsnFrame = False snres@lbLabelBarOn = False ; turn off label bar snres@cnInfoLabelOn = False ; turn off info label snres@cnLineLabelsOn = False ; turn off contour labels snres@cnLinesOn = False ; no contour lines snres@cnConstFEnableFill= True ; allow constant values (1=sig) to fill snres@cnFillOn = True ; turn on fill ;snres@cnMonoFillPattern = True ; single fill pattern snres@cnFillColor = "red" ; snres@cnFillPattern = 17 ; stipling snres@cnFillDotSizeF = 0.003 snres@cnFillScaleF = 0.5 snres@cnFillMode = "areafill" ;snres@cnFillMode = "cellfill" ; try this to see the true extent of the s/n snres@cnMonoFillColor = True snres@cnConstFLabelOn = False ;---------------------------------------------------------------------- ; Loop through six PRISM months and 12 WRF files to calculate difference ;---------------------------------------------------------------------- plots = new(6,graphic) labels = new(6,string) do nfp=0,nfilp-1 yyyymm = toint( str_get_field(filp(nfp),5,"_.") ) ; parse PRISM name for time yyyy = yyyymm/100 ; current year as integer mm = yyyymm-(yyyy*100) ; current month labels(nfp) = yyyy+sprinti("%0.2i",mm) fp = addfile(dirp+filp(nfp), "r") XP = fp->z ; original PRISM variable (all 'z') xp = ESMF_regrid_with_weights(XP,wgtFilePath,False); create new variable on WRF grid printVarSummary(xp) ;----Loop through 12 WRF ensembles Wrfdiff = new((/12,192,132/),float,1e20) do nfc=0,ndirPD-1 ; matching WRF file filw = "avg_"+varw+"_"+yyyy+"-"+sprinti("%0.2i",mm)+".nc" ; Present-Day average temperature, snowc, snow print(filw) ;filw = "total_"+varw+"_"+yyyy+"-"+sprinti("%0.2i",mm)+".nc" ; total accumulated precip if (isfilepresent(dirPD+dirPDs(nfc)+"/nc/"+filw)) then fPD = addfile(dirPD+dirPDs(nfc)+"/nc/"+filw, "r") ; Open WRF Present-Day file (if present) xPD = fPD->$varw$ ; WRF present day monthly variable fDF = addfile(dirDF+dirDFs(nfc)+"/nc/"+filw, "r") xDF = fDF->$varw$ printVarSummary(xPD) ;---------------------------------------------------------------------- ; Subtract PRISM from WRF, masked and average by LC type conversion ;---------------------------------------------------------------------- ;---- Create new variable to hold differences between WRF and Prism xdiff = xPD(0,:,:) xdiff_MF2GR = xPD(0,:,:) xdiff_EN2GR = xPD(0,:,:) xdiff_DB2GR = xPD(0,:,:) xdiff_UR2GR = xPD(0,:,:) xdiff_CW2GR = xPD(0,:,:) xdiff = xPD(0,:,:)-xp ; Avg monthly WRF minus Prism for month i xdiff_MF2GR = mask(xdiff,ludiff,10) ; Bias of Mixed Forest to Grassland (15->5, 10) xdiff_EN2GR = mask(xdiff,ludiff,9) ; Bias of Evergreen Needleleaf to grass (14->5, 9) xdiff_DB2GR = mask(xdiff,ludiff,6) ; Bias of Decid. Broadleaf to Grassland (11->5, 6) xdiff_UR2GR = mask(xdiff,ludiff,-4) ; bias of Urban to Grassland (1 ->5, -4) xdiff_CW2GR = mask(xdiff,ludiff,1) ; bias of Crop/Wood to Grassland (6 ->5, 1) xdiff_CRGR = mask(xdiff,luPD,5) ; bias of Cropland/Grassland Mosaic (5) ;---- Create new variables to hold biases by LC type conversion bias_CRGR = avg(xdiff_CRGR) bias_MF2GR = avg(xdiff_MF2GR)-bias_CRGR bias_EN2GR = avg(xdiff_EN2GR)-bias_CRGR bias_DB2GR = avg(xdiff_DB2GR)-bias_CRGR bias_UR2GR = avg(xdiff_UR2GR)-bias_CRGR bias_CW2GR = avg(xdiff_CW2GR)-bias_CRGR print("bias_CRGR= "+bias_CRGR) print("bias_MF2GR= "+bias_MF2GR) print("bias_EN2GR= "+bias_EN2GR) print("bias_DB2GR= "+bias_DB2GR) print("bias_UR2GR= "+bias_UR2GR) print("bias_CW2GR= "+bias_CW2GR) ;---------------------------------------------------------------------- ; Apply PRISM biases to Present-Day minus Deforested for each LC ; Conversion type ;---------------------------------------------------------------------- LCdiff = xPD(0,:,:) printVarSummary(LCdiff) printVarSummary(xDF) LCdiff = xPD(0,:,:)-xDF(0,:,:) LCdiff = where(ludiff.eq.10,LCdiff-bias_MF2GR,LCdiff) LCdiff = where(ludiff.eq.9,LCdiff-bias_EN2GR,LCdiff) LCdiff = where(ludiff.eq.6,LCdiff-bias_DB2GR,LCdiff) LCdiff = where(ludiff.eq.-4,LCdiff-bias_UR2GR,LCdiff) LCdiff = where(ludiff.eq.1,LCdiff-bias_CW2GR,LCdiff) LCdiff = mask(LCdiff,lm,1) LCdiff@lat2d = wrfPD->XLAT(0,:,:) LCdiff@lon2d = wrfPD->XLONG(0,:,:) Wrfdiff(nfc,:,:)=LCdiff ; assign xdiff to WrfPrdiff variable end if ; isfilepresent end do ; nfc Wrfdiff!0 = "ensemble" ;---------------------------------------------------------------------- ; Average bias corrected Present-Day minus Deforested for region and ; for each land cover type conversion. Plot the regional average ;---------------------------------------------------------------------- Wrfdiff_mean = dim_avg_n_Wrap(Wrfdiff,0) ; ensemble average printVarSummary(Wrfdiff_mean) print("min Wrfdiff_mean= "+min(Wrfdiff)) print("max Wrfdiff_mean= "+max(Wrfdiff)) Wrfdiff_std = dim_stddev_n_Wrap(Wrfdiff,0) ; ensemble standard deviation printVarSummary(Wrfdiff_std) print("min Wrfdiff_std= "+min(Wrfdiff_std)) print("max Wrfdiff_std= "+max(Wrfdiff_std)) Wrfdiff_sn = Wrfdiff_mean ; copy metadata to signal/noise (sn) Wrfdiff_sn = Wrfdiff_mean/Wrfdiff_std ; calculate signal to noise printVarSummary(Wrfdiff_sn) print("min SN= "+min(Wrfdiff_sn)) print("max SN= "+max(Wrfdiff_sn)) Wrfdiff_sn2 = Wrfdiff_sn ; copy metadata to sn >2 |<-2 Wrfdiff_sn2 = where(abs(Wrfdiff_sn).lt.2,1,Wrfdiff@_FillValue) Wrfdiff_sn2@lat2d = wrfPD->XLAT(0,:,:) Wrfdiff_sn2@lon2d = wrfPD->XLONG(0,:,:) printVarSummary(Wrfdiff_sn2) print("min Wrfdiff_sn2= "+min(Wrfdiff_sn2)) print("max Wrfdiff_sn2= "+max(Wrfdiff_sn2)) plot_diff = gsn_csm_contour_map(wks,Wrfdiff_mean,res) plot_sn = gsn_csm_contour(wks,Wrfdiff_sn2,snres) overlay(plot_diff,plot_sn) plots(nfp)=plot_diff end do ; nfp ;---------------------------------------------------------------------- ; Panel Plotting options section ;---------------------------------------------------------------------- pres = True pres@gsnMaximize = True pres@gsnPanelLabelBar = True pres@gsnPanelYWhiteSpace= 5 pres@gsnPanelXWhiteSpace= 5 pres@txString = "Present-Day minus Deforested" pres@gsnPanelFigureStrings = labels pres@gsnPanelFigureStringsFontHeightF = 0.006 pres@gsnPanelFigureStringsPerimOn = True gsn_panel(wks,plots,(/1,6/),pres) do i = 0, 5 draw(plots(i)) frame(wks) end do end