;plots 3D(time,ns,we) field difference for list of years ;and months. no vector fields. ;for variable specified. plots current, future, absolute diff, percent diff ;DOES NOT test significance of differences ;does average of the field only for years and months specified ;uses monthly average data! ;for native gridded data only ;*************************** ;NOTE!... ;will do all RCM+GCM combinations found in the field's folder ;for for the specified resolution and RCP ;data must be organized such that all simulations for a given field are ;in their own folder ;*************************** ;MSB 2017 load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_csm.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl" load "/glade/u/home/bukovsky/ncl/cordex/functions_cordex/cordex_read_data_funct.ncl" load "/glade/u/home/bukovsky/ncl/cordex/functions_cordex/cordex_plotting.ncl" load "/glade/u/home/bukovsky/ncl/cordex/functions_cordex/cordex_convert.ncl" load "/glade/u/home/bukovsky/ncl/cordex/functions_cordex/cordex_stats.ncl" load "/glade/u/home/bukovsky/ncl/cordex/functions_cordex/cordex_bootstrapping.ncl" ;load "/glade/u/home/bukovsky/ncl/phd/my_ncl_funct/calc_climate.ncl" ;load "/glade/u/home/bukovsky/ncl/phd/my_ncl_funct/calc_stat.ncl" ;load "/glade/u/home/bukovsky/ncl/phd/my_ncl_funct/my_plot_functions.ncl" ;load "/glade/u/home/bukovsky/ncl/phd/my_ncl_funct/convert.ncl" ;load "/glade/u/home/bukovsky/ncl/phd/WRF/my_wrf_functs/convertWRF.ncl" ;load "/glade/u/home/bukovsky/ncl/phd/WRF/my_wrf_functs/calc_wrf_msb.ncl" ;load "/glade/u/home/bukovsky/ncl/my_narccap_functs/narccap_plotting.ncl" ;load "/glade/u/home/bukovsky/ncl/my_narccap_functs/narccap_functs.ncl" begin ;---------------------------------------------------------------------- futures = (/"rcp45","rcp85"/) resolutions = (/"22","44","11"/) ;---------------------------------------------------------------------- ;dir = "/data/" dir = "/esgf/" future = futures(1) ;resolution = resolutions(1) resolution = resolutions(0) ; number of years between 1 and 2 must be the same for bootstrapping function!! ; must have more than 1 year too yrstart1 = "1951" yrend1 = "1999" yrstart2 = "2051" yrend2 = "2099" ;MONTHS AS INTEGERS FOR AVERAGING ;must be the same ;mons = (/12,1,2/) ;mons = (/3,4,5/) mons = (/6,7,8/) ;mons = (/9,10,11/) ;mons = (/1,2,3,4,5,6,7,8,9,10,11,12/) ;VARIABLE field = "pr" ;field = "tas" ;REGION, IF ONLY PART OF DOMAIN WANTED, SPECIFY AS TRUE AND GIVE COORDS ;region = False region = True ;CORDEX COMMON minlat = 18. maxlat = 60. minlon = -165. maxlon = -65. ;NARCCAP COMMON ;minlat = 22.75 ;maxlat = 69.75 ;minlon = -149.25 ;maxlon = -43.75 ;FACETS ;minlat = 12.5 ;maxlat = 50. ;minlon = -123.7 ;maxlon = -47.6 ;CUS ;minlat = 32. ;maxlat = 49. ;minlon = -105. ;maxlon = -85. ;CUS+ ;minlat = 30. ;maxlat = 50. ;minlon = -105. ;maxlon = -85. ;NAM ;minlat = 20. ;maxlat = 40. ;minlon = -116. ;maxlon = -104. ;NAM2 ;minlat = 20. ;maxlat = 39. ;minlon = -115. ;maxlon = -103. ;NAM3 ;minlat = 19. ;maxlat = 42. ;minlon = -117. ;maxlon = -101. ;NE CANADA ;minlat = 50. ;maxlat = 75. ;minlon = -110. ;maxlon = -40. ;10km Western Domain ;minlat = 21.0 ;maxlat = 55.0 ;minlon = -124.0 ;maxlon = -85.0 ;Mid Atlantic: WRFG,RCM3,MM5I,HRM3 ;minlat = 34.5 ;maxlat = 40. ;minlon = -81. ;maxlon = -73.5 ;Mid Atlantic: CRCM, ECP2 ;minlat = 34.5 ;maxlat = 39.5 ;minlon = -81. ;maxlon = -73. ;Mid Atlantic: TMSL ;minlat = 34.0 ;maxlat = 40.5 ;minlon = -80. ;maxlon = -74.5 ;South CUS ;minlat = 32. ;maxlat = 40. ;minlon = -105. ;maxlon = -91. ;East South Central ;minlat = 20. ;maxlat = 53. ;minlon = -106. ;maxlon = -65. ;RASTER MAP PLOTS? ;rasterize = True rasterize = False ;---------------------------------------------------------------------- ;---------------------------------------------------------------------- ;A FEW NEEDED THINGS years1 = ispan(stringtoint(yrstart1),stringtoint(yrend1),1) years2 = ispan(stringtoint(yrstart2),stringtoint(yrend2),1) smons = sprinti("%0.2i",mons) ;for plot titles ;for output plot file name mletters = (/"J","F","M","A","M","J","J","A","S","O","N","D"/) nmons = dimsizes(mons) if(nmons.eq.12)then fmons = "ANN" else do i = 0,nmons-1 if(i.eq.0)then fmons = mletters(mons(i)-1) else fmons = fmons+mletters(mons(i)-1) end if end do end if if(region.eq.False) minlat = -999. maxlat = -999. minlon = -999. maxlon = -999. end if ;---------------------------------------------------------------------- ;---------------------------------------------------------------------- ;FIND FILES TO PLOT ;figure out what models we will be pulling and plotting ;file name structure: field.rcp.GCM.RCM.mon.resolution.raw.nc if(dir.eq."/data/") fhist = systemfunc("ls "+dir+field+"/"+field+".hist.*.*.mon.*-"+resolution+".raw.nc") ffuture = systemfunc("ls "+dir+field+"/"+field+"."+future+".*.*.mon.*-"+resolution+".raw.nc") end if if(dir.eq."/esgf/") fhist = systemfunc("ls "+dir+field+"/"+field+"_*-"+resolution+"_*_historical_*_mon_*.nc") ffuture = systemfunc("ls "+dir+field+"/"+field+"_*-"+resolution+"_*_"+future+"_*_mon_*.nc") end if if(.not.isdefined("fhist"))then print("NEED TO DEFINE NAMING STRUCTURE FOR THIS DIRECTORY OF DATA") exit end if nf0 = dimsizes(fhist) nf1 = dimsizes(ffuture) print(fhist) print(ffuture) print(" ") print("Find corresponding historical and future simulations.") print(" ") rcmGCM0 = new(nf0,string) rcmGCM1 = new(nf1,string) if(dir.eq."/data/") do i = 0,nf0-1 str0 = str_split(fhist(i),".") rcmGCM0(i) = str0(2)+"."+str0(3) end do do i = 0,nf1-1 str1 = str_split(ffuture(i),".") rcmGCM1(i) = str1(2)+"."+str1(3) end do end if if(dir.eq."/esgf/") do i = 0,nf0-1 str0 := str_split(fhist(i),"_") rcmGCM0(i) = str0(5)+"."+str0(2) end do do i = 0,nf1-1 str1 := str_split(ffuture(i),"_") rcmGCM1(i) = str1(5)+"."+str1(2) end do end if allrcmgcm = array_append_record(rcmGCM0,rcmGCM1,0) allunique = get_unique_values(allrcmgcm) ; print(allrcmgcm) ; print(allunique) allunique@_FillValue = default_fillvalue(typeof(allunique)) nu = dimsizes(allunique) fhtmp = new( (/nu,20/),string) fftmp = new( (/nu,20/),string) do i = 0,nu-1 iwantc = ind(rcmGCM0.eq.allunique(i)) iwantf = ind(rcmGCM1.eq.allunique(i)) ; print(i) ; print(iwantc) ; print(iwantf) ; print(.not.ismissing(iwantc)) ; print(.not.all(ismissing(iwantf))) if((all(ismissing(iwantc))).or.(all(ismissing(iwantf))))then allunique(i) = allunique@_FillValue else ; print(iwantc+" "+fhist(iwantc)) ; print(iwantf+" "+ffuture(iwantf)) nh = dimsizes(iwantc) nf = dimsizes(iwantf) fhtmp(i,0:nh-1) = fhist(iwantc) fftmp(i,0:nf-1) = ffuture(iwantf) end if delete(iwantc) delete(iwantf) end do inmu = ind(.not.ismissing(allunique)) delete(ffuture) delete(fhist) fhist = fhtmp(inmu,:) ffuture = fftmp(inmu,:) delete(fhtmp) delete(fftmp) ; print(fhist) ; print(ffuture) delete(str0) delete(str1) delete(rcmGCM0) delete(rcmGCM1) delete(nf0) delete(nf1) nsim = dimsizes(allunique(inmu)) ;get names of RCM-GCM combinations rcm_gcm = new(nsim,string) rcm_gcm = allunique(inmu) rcm_gcm = str_sub_str(rcm_gcm,".","+") rcm_gcm = str_sub_str(rcm_gcm,"_","+") print(rcm_gcm) ;---------------------------------------------------------------------- ;---------------------------------------------------------------------- ;SET PLOTTING STUFF THAT DOES NOT NEED TO BE IN THE LOOP ;OPEN OUTPUT DEVICE ;pltype = "pdf" pltype = "png" ;plname = "cordex_"+nsim+"sim_"+future+"_"+resolution+"_"+yrstart1+"-"+yrend1+"_"+yrstart2+"-"+yrend2+"_"+fmons+"_"+field+"_nosig" ;wks = gsn_open_wks(pltype,plname) ;define 2 colormaps, one for current/future, one for differences ;colormap1 = read_colormap_file("BlAqGrYeOrReVi200") ;colormap2 = read_colormap_file("testcmap") colormap1 = read_colormap_file("GMT_wysiwygcont") colormap2 = read_colormap_file("NCV_jaisnd") ;---------------------------------------------------------------------- ;PLOT ;GENERAL PLOTTING STUFF opts = True opts@mpFillOn = False opts@mpOutlineDrawOrder = "PostDraw" opts@mpOutlineBoundarySets = "GeophysicalAndUSStates" opts@mpPerimDrawOrder = "PostDraw" opts@gsnAddCyclic = False ;for regional opts@cnSpanFillPalette = True opts@gsnMaximize = True opts@cnFillOn = True opts@cnLineLabelsOn = False opts@cnLinesOn = False if(rasterize)then opts@cnFillMode = "CellFill" end if opts@lbLabelBarOn = True ;common label bar opts@lbBoxMinorExtentF = 0.15 opts@lbLabelAngleF = -90. opts@lbLabelAutoStride = True opts@lbLabelFontHeightF = .014 opts@lbTitleOn = True opts@lbTitleFontHeightF = .015 opts@lbTitlePosition = "Bottom" opts@lbTitleDirection = "Across" opts@pmLabelBarOrthogonalPosF = -0.02 ;turn on/off lat/lon labeling opts@pmTickMarkDisplayMode = "Never" ;opts@gsnDraw = False ;opts@gsnFrame = False ;---------------------------------------------------------------------- ;RESOURCES FOR STIPPLING ;ores = True ;ores = set_stipple_opts(ores,"LT","white") ;---------------------------------------------------------------------- ;---------------------------------------------------------------------- ;---------------------------------------------------------------------- ;---------------------------------------------------------------------- do n = 0,nsim-1 ;big loop across simulations print(" ") print(" ") print("DOING: "+rcm_gcm(n)) print(" ") print(" ") ;---------------------------------------------------------------------- ;---------------------------------------------------------------------- ;GET VARIABLE FROM FILE inmh = ind(.not.ismissing(fhist(n,:))) inmf = ind(.not.ismissing(ffuture(n,:))) ; print(inmh) ; print(fhist(n,:)) ; print(inmf) ; print(ffuture(n,:)) var1 = cordex_get_data_fromfile(fhist(n,inmh),yrstart1,yrend1,mons,field, \ minlat,maxlat,minlon,maxlon) var2 = cordex_get_data_fromfile(ffuture(n,inmh),yrstart2,yrend2,mons,field, \ minlat,maxlat,minlon,maxlon) delete(inmh) delete(inmf) date1 = tointeger(var1&time);comes off of variable as not an integer, needs to be date2 = tointeger(var2&time) lat = var1@lat2d lon = var1@lon2d ;printVarInfo(date1,"date1") ;printVarInfo(date2,"date2") ;printVarInfo(lat,"lat") ;printVarInfo(lon,"lon") ;printVarInfo(var1,"var1") ;printVarInfo(var2,"var2") ;---------------------------------------------------------------------- dims1 = dimsizes(var1) dims2 = dimsizes(var2) nt1 = dims1(0) nt2 = dims2(0) nlat = dims1(1) nlon = dims1(2) if(field.eq."pr".or.field.eq."prc".or.field.eq."evps")then var1 = var1 > 0.0 ;make sure all values are gt zero var2 = var2 > 0.0 ; change rate from mm/s to mm/day var1 = var1*3600*24 var1@units = "mm/day" var2 = var2*3600*24 var2@units = "mm/day" end if if(field.eq."tas")then var1 = var1-273.15 var2 = var2-273.15 var1@units = "C" var2@units = "C" end if if(field.eq."huss")then var1 = var1 > 0.0 var2 = var2 > 0.0 var1 = var1*1000. var2 = var2*1000. var1@units = "g/kg" var2@units = "g/kg" end if ;printVarInfo(var1,"var1") ;printVarInfo(var2,"var2") ;---------------------------------------------------------------------- ;---------------------------------------------------------------------- ;CALCULATE STUFF, THEIR DIFFERENCE AND PERCENT GROWTH ;AVERAGE varavg1 = dim_avg_n_Wrap(var1,0) varavg2 = dim_avg_n_Wrap(var2,0) ;printVarInfo(varavg1,"varavg1") ;printVarInfo(varavg2,"varavg2") davg = varavg2 davg = 0.0 davg = varavg2-varavg1 pgavg = percent_growth(varavg1,varavg2) avg1 = avg(var1) avg2 = avg(var2) ;printVarInfo(davg,"davg") ;printVarInfo(pgavg,"pgavg") print("var1 avg: "+avg1) print("var2 avg: "+avg2) print("davg avg: "+avg(davg)) print("pgavg avg: "+avg(pgavg)) ;---------------------------------------------------------------------- ;---------------------------------------------------------------------- ;CREATE PLOTS! ;set in projection functions ;opts@tfDoNDCOverlay = True ;only set to do if plotting with native projection ;opts@tfDoNDCOverlay = False res = opts res = set_forced_proj_2d(var1,opts,"Satellite") ;rese= set_map_proj(var1,opts) ;sets native projections ;---------------------------------------------------------------------- ;PLOT AVERAGE CURRENT res@cnFillPalette = colormap1 ;set color map res@cnLevelSelectionMode = "AutomaticLevels" if(field.eq."pr".or.field.eq."prc".or.field.eq."prec")then res@cnLevelSelectionMode = "ManualLevels" res@cnMinLevelValF = 0.0 ; res@cnMaxLevelValF = 6.0 ; res@cnLevelSpacingF = 0.2 res@cnMaxLevelValF = 10.0 res@cnLevelSpacingF = 0.5 res@cnFillPalette = colormap1(::-1,:) end if if(field.eq."tas".or.field.eq."tmin".or.field.eq."tmax")then res@cnLevelSelectionMode = "ManualLevels" res@cnMinLevelValF = -40. res@cnMaxLevelValF = 40. res@cnLevelSpacingF = 5. ; res@cnMinLevelValF = 0. ; res@cnMaxLevelValF = 40. ; res@cnLevelSpacingF = 2. end if if(field.eq."huss")then res@cnLevelSelectionMode = "ManualLevels" res@cnMinLevelValF = 0. res@cnMaxLevelValF = 20. res@cnLevelSpacingF = 1. res@cnFillPalette = colormap1(::-1,:) end if res@lbTitleString = "("+var1@units+")" res@tiMainString = "~Z50~"+var1@long_name+" "+rcm_gcm(n)+" 0."+resolution+" deg" res@gsnRightString = "Avg = "+avg(varavg1) res@gsnLeftString = fmons+" "+yrstart1+"-"+yrend1 plname = rcm_gcm(n)+"_"+resolution+"_"+yrstart1+"-"+yrend1+"_"+fmons+"_"+field+"_hist" wks = gsn_open_wks(pltype, plname) plot = gsn_csm_contour_map(wks,varavg1,res) ;---------------------------------------------------------------------- ;PLOT AVERAGE FUTURE getvalues plot@contour "cnLevels" : plevels end getvalues ;print(plevels) delete(res@cnLevelSelectionMode) res@cnLevelSelectionMode = "ExplicitLevels" res@cnLevels = plevels res@gsnLeftString = fmons+" "+yrstart2+"-"+yrend2 res@gsnRightString = "Avg = "+avg(varavg2) plname = rcm_gcm(n)+"_"+resolution+"_"+yrstart2+"-"+yrend2+"_"+fmons+"_"+field+"_future" wks = gsn_open_wks(pltype, plname) plot = gsn_csm_contour_map(wks,varavg2,res) ;---------------------------------------------------------------------- ;ABSOLUTE DIFFERENCE delete(res@cnLevelSelectionMode) delete(res@cnLevels) res@gsnLeftString = fmons+" "+yrstart2+"-"+yrend2+" - "+yrstart1+"-"+yrend1 res@gsnRightString = "Avg = "+avg(davg) delete(res@cnFillPalette) if(field.eq."pr".or.field.eq."prc".or.field.eq."huss")then res@cnFillPalette = colormap2(::-1,:) ;set color map reversed end if if(field.eq."tas".or.field.eq."tmin".or.field.eq."tmax")then res@cnFillPalette = colormap2 end if if(.not.isatt(res,"cnFillPalette"))then res@cnFillPalette = colormap2 end if symMinMaxPlt(davg,20,False,res) if(field.eq."pr".or.field.eq."prc")then res@cnLevelSelectionMode = "ManualLevels" res@cnMinLevelValF = -2.6 res@cnMaxLevelValF = 2.6 res@cnLevelSpacingF = 0.2 end if if(field.eq."huss")then res@cnLevelSelectionMode = "ManualLevels" res@cnMinLevelValF = -5. res@cnMaxLevelValF = 5. res@cnLevelSpacingF = 0.2 if(any(mons.eq.12))then res@cnMinLevelValF = -2. res@cnMaxLevelValF = 2. res@cnLevelSpacingF = 0.1 end if end if if(field.eq."tas".or.field.eq."tmin".or.field.eq."tmax")then res@cnLevelSelectionMode = "ExplicitLevels" ; res@cnLevels = (/-3,-2,-1,0,0.5,1,1.5,2,3,4,5,7,10/) res@cnLevels = (/-3,-2,-1,0,1,2,3,4,5,7,9,12/) lminus = (/2,40,100,125/) ; lplus = round(fspan(130,255,10),3) lplus = round(fspan(130,255,9),3) res@cnFillColors = array_append_record(lminus,lplus,0) end if plname = rcm_gcm(n)+"_"+resolution+"_"+yrstart2+"-"+yrend2+"_"+fmons+"_"+field+"_absdiff" wks = gsn_open_wks(pltype, plname) plot = gsn_csm_contour_map(wks,davg,res) if(field.eq."tas")then delete(res@cnLevels) delete(res@cnFillColors) end if ;---------------------------------------------------------------------- ;PERCENT DIFFERENCE res@gsnLeftString = fmons+" "+yrstart1+"-"+yrend1+" vs. "+yrstart2+"-"+yrend2 res@gsnRightString = "Avg = "+avg(pgavg) res@lbTitleString = "(%)" symMinMaxPlt(pgavg,20,False,res) if(field.eq."pr".or.field.eq."prc".or.field.eq."prec")then res@cnLevelSelectionMode = "ManualLevels" res@cnMinLevelValF = -100. res@cnMaxLevelValF = 100. res@cnLevelSpacingF = 10. ; res@cnMinLevelValF = -50. ; res@cnMaxLevelValF = 50. ; res@cnLevelSpacingF = 5. end if if(field.eq."tas".or.field.eq."tmin".or.field.eq."tmax")then res@cnLevelSelectionMode = "ManualLevels" res@cnMinLevelValF = -200. res@cnMaxLevelValF = 200. res@cnLevelSpacingF = 20. end if if(field.eq."huss")then res@cnLevelSelectionMode = "ManualLevels" res@cnMinLevelValF = -30. res@cnMaxLevelValF = 30. res@cnLevelSpacingF = 2. if(any(mons.eq.12))then res@cnMinLevelValF = -50. res@cnMaxLevelValF = 50. res@cnLevelSpacingF = 5. end if end if plname = rcm_gcm(n)+"_"+resolution+"_"+yrstart2+"-"+yrend2+"_"+fmons+"_"+field+"_pctdiff" wks = gsn_open_wks(pltype, plname) plot = gsn_csm_contour_map(wks,pgavg,res) ;---------------------------------------------------------------------- ;---------------------------------------------------------------------- ;---------------------------------------------------------------------- ;END DO LOOP ACROSS SIMULATIONS units = var1@units delete(var1) delete(var2) delete(date1) delete(date2) delete(lat) delete(lon) delete(varavg1) delete(varavg2) delete(davg) delete(pgavg) delete(res) end do ;---------------------------------------------------------------------- ;---------------------------------------------------------------------- print("DONE") end