; Julia 22/01/2014: This script calculates the following precipitation indices; R95pTOT, R99pTOT, PRCPTOT. We are using the historical runs as the baseline period. 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 "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" load "/scratch/y98/shared/ncl-functions/cordex_grid_plot_functions.ncl" load "/scratch/y98/shared/ncl-functions/max_run_length.ncl" begin ;---------------INPUTS------------------------------------- x = 100 y = 100 R95pTOT_era = new((/29,x,y/),float) R95pTOT_bom = new((/29,x,y/),float) R99pTOT_era = new((/29,x,y/),float) R99pTOT_bom = new((/29,x,y/),float) PRCPTOT_era = new((/29,x,y/),float) PRCPTOT_bom = new((/29,x,y/),float) SPII_bom = new((/29,x,y/),float) SPII_era = new((/29,x,y/),float) RX5d_era = new((/29,x,y/),float) RX5d_bom = new((/29,x,y/),float) year = (/1981,1982,1983,1984,1985,1986,1987,1988,1989,1990,1991,1992,1993,1994,1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007,2008,2009/) name = "ERA" file_suffix = "tzadjust" domain = "d03" filo = name + "_" + domain + "_indices_prec.nc" diro = "Processed/" diro_bom = "/group/y98/julia/BOM/Processed/" ;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 delete(f) ;load all the BOM files and work out the 95th and 99th percentiles. We just need to keep the percentiles right now files_bom = systemfunc("ls " + diro_bom + "BOM_" + domain + "_*.nc") f_bom = addfiles(files_bom + ".nc", "r") ListSetType (f_bom, "cat") PPT_bom = f_bom[:]->DLY_PPT(:,:,:) delete(f_bom) dsizes_array = dimsizes(PPT_bom) PPT_bom_nr = PPT_bom(:,10:(dsizes_array(1)-11),10:(dsizes_array(2)-11)) delete(PPT_bom) delete(dsizes_array) ppt_95_bom_baseline = new((/1,x,y/),float) ppt_99_bom_baseline = new((/1,x,y/),float) counteri = 0 do i = 0,x-1,1 counterj = 0 do j = 0,y-1,1 ;rank precip values ppt_bom = PPT_bom_nr(:,counteri,counterj) qsort(ppt_bom) ;calculate N, where N is the number of days with ppt ge 1 N = num(ppt_bom.ge.1) Ninv = num(ppt_bom.lt.1) ;calculate percentile rank rank95_ppt_bom = 0.95 * N rank99_ppt_bom = 0.99 * N ; we need to round the percentile rank to a whole number and then we add Ninv, to capture the correct index value accounting for no rain days. rank95_ppt_bom_round = round(rank95_ppt_bom,3) + Ninv rank99_ppt_bom_round = round(rank99_ppt_bom,3) + Ninv delete(rank95_ppt_bom) delete(rank99_ppt_bom) ;now we are using these integers to extract the percentile value ppt95_val = ppt_bom(rank95_ppt_bom_round-1) ppt99_val = ppt_bom(rank99_ppt_bom_round-1) delete(rank95_ppt_bom_round) delete(rank99_ppt_bom_round) ;and finally we return these values to the main array ppt_95_bom_baseline(:,counteri,counterj) = ppt95_val ppt_99_bom_baseline(:,counteri,counterj) = ppt99_val delete(ppt95_val) delete(ppt99_val) counterj = counterj + 1 end do counteri = counteri + 1 end do ;print(ppt_95_bom_baseline) delete(PPT_bom_nr) do k = 0,28,1 files_input_era = systemfunc("ls " + diro + name + "_" + domain + "_" + year(k) + "_tzadjust.nc") files_input_bom = systemfunc("ls " + diro_bom + "BOM_" + domain + "_" +year(k) + ".nc") ;------load files, a year at a time and work out our indices --------------- f_era = addfiles(files_input_era + ".nc", "r") ListSetType (f_era, "join") f_bom = addfiles(files_input_bom + ".nc", "r") ListSetType (f_bom, "join") ; -------- now extract the vairables we need PPT_era = f_era[:]->DLY_PPT PPT_bom = f_bom[:]->DLY_PPT delete(f_bom) delete(f_era) ;mask water and remove relaxation zone dsizes_array = dimsizes(landmask) 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) ; we need to remove the relaxation zone from the BOM data too PPT_bom_nr = PPT_bom(:,10:(dsizes_array(0)-11),10:(dsizes_array(1)-11)) ppt_era_mask = mask(PPT_era,landmask_nr,1) ppt_bom_mask = mask(PPT_bom_nr,landmask_nr,1) delete(PPT_era) delete(PPT_bom) delete(PPT_bom_nr) ; we are now going to work out the 5 days running mean. ppt_era_runave = runave_n(ppt_era_mask,5,0,0) ppt_bom_runave = runave_n(ppt_bom_mask,5,0,0) ppt_era_runave_max = dim_max_n(ppt_era_runave,0) ppt_bom_runave_max = dim_max_n(ppt_bom_runave,0) delete(ppt_era_runave) delete(ppt_bom_runave) RX5d_era(k,:,:) = ppt_era_runave_max * 5 RX5d_bom(k,:,:) = ppt_bom_runave_max * 5 ;we are now working out the yearly SPII index, and PRCPTOT (which is annual count of wet days) ppt_era_adj = where(ppt_era_mask.lt.1,0,ppt_era_mask) RR_era = dim_sum_n(ppt_era_adj,0) W_era = dim_num_n(ppt_era_adj.ge.1,0) ;these extra steps are to remove the zeros so we can divide by the variable W_era_fv = where(W_era.ne.0,W_era,W_era@FillValue) spii_era = RR_era/W_era_fv ;spii_era = where(ismissing(spii_era),0,y) delete(ppt_era_adj) ppt_bom_adj = where(ppt_bom_mask.lt.1,0,ppt_bom_mask) RR_bom = dim_sum_n(ppt_bom_adj,0) W_bom = dim_num_n(ppt_bom_adj.ge.1,0) W_bom_fv = where(W_bom.ne.0,W_bom,W_bom@FillValue) spii_bom = RR_bom/W_bom_fv ;spii_bom = where(ismissing(spii_bom),0,y) delete(ppt_bom_adj) SPII_era(k,:,:) = spii_era SPII_bom(k,:,:) = spii_bom PRCPTOT_era(k,:,:) = W_era PRCPTOT_bom(k,:,:) = W_bom delete(spii_era) delete(spii_bom) delete(W_era) delete(W_bom) ;we are now using some of the above variables to work out R95pTOT and R99pTOT counterl = 0 do l = 0,x-1,1 counterm = 0 do m = 0,y-1,1 ;extract the l,m percentile values ppt_95th = ppt_95_bom_baseline(:,l,m) ppt_99th = ppt_99_bom_baseline(:,l,m) ;now use these values to zero values below threshold ppt_era_95 = where(ppt_era_mask(:,counterl,counterm).lt.ppt_95_bom_baseline(:,counterl,counterm),0,ppt_era_mask(:,counterl,counterm)) ppt_era_99 = where(ppt_era_mask(:,counterl,counterm).lt.ppt_99_bom_baseline(:,counterl,counterm),0,ppt_era_mask(:,counterl,counterm)) ppt_bom_95 = where(ppt_bom_mask(:,counterl,counterm).lt.ppt_95_bom_baseline(:,counterl,counterm),0,ppt_bom_mask(:,counterl,counterm)) ppt_bom_99 = where(ppt_bom_mask(:,counterl,counterm).lt.ppt_99_bom_baseline(:,counterl,counterm),0,ppt_bom_mask(:,counterl,counterm)) ;we then sum all the precipitation above the thresholds ppt_era_95_sum = dim_sum_n(ppt_era_95,0) ppt_era_99_sum = dim_sum_n(ppt_era_99,0) ppt_bom_95_sum = dim_sum_n(ppt_bom_95,0) ppt_bom_99_sum = dim_sum_n(ppt_bom_99,0) delete(ppt_era_95) delete(ppt_era_99) delete(ppt_bom_95) delete(ppt_bom_99) ;and return the values to the array R95pTOT_era(k,counterl,counterm) = ppt_era_95_sum R99pTOT_era(k,counterl,counterm) = ppt_era_99_sum R95pTOT_bom(k,counterl,counterm) = ppt_bom_95_sum R99pTOT_bom(k,counterl,counterm) = ppt_bom_99_sum delete(ppt_era_95_sum) delete(ppt_era_99_sum) delete(ppt_bom_95_sum) delete(ppt_bom_99_sum) counterm = counterm + 1 end do counterl = counterl + 1 end do delete(ppt_era_mask) delete(ppt_bom_mask) end do ; now that we have our annual data, we need to take the clim means R95pTOT_era_mean = dim_avg_n(R95pTOT_era,0) R95pTOT_bom_mean = dim_avg_n(R95pTOT_bom,0) R99pTOT_era_mean = dim_avg_n(R99pTOT_era,0) R99pTOT_bom_mean = dim_avg_n(R99pTOT_bom,0) PRCPTOT_era_mean = dim_avg_n(PRCPTOT_era,0) PRCPTOT_bom_mean = dim_avg_n(PRCPTOT_bom,0) SPII_era_mean = dim_avg_n(SPII_era,0) SPII_bom_mean = dim_avg_n(SPII_bom,0) RX5d_era_mean = dim_avg_n(RX5d_era,0) RX5d_bom_mean = dim_avg_n(RX5d_bom,0) ;we need to mask PRCPTOT again PRCPTOT_bom_mean_mask = mask(PRCPTOT_bom_mean,landmask_nr,1) PRCPTOT_era_mean_mask = mask(PRCPTOT_era_mean,landmask_nr,1) ;---- write to a netcdf file for further processing --------- system("/bin/rm -f " + filo) fout = addfile(diro+filo,"c") fout->lat = (/lat/) fout->lon = (/lon/) fout->SPII_era = (/SPII_era_mean/) fout->SPII_bom = (/SPII_bom_mean/) fout->RX5d_era = (/RX5d_era_mean/) fout->RX5d_bom = (/RX5d_bom_mean/) fout->R95pTOT_era = (/R95pTOT_era_mean/) fout->R95pTOT_bom = (/R95pTOT_bom_mean/) ;---- now we plot ;fi = new((/3,x,y/),float) ; ;fi(0,:,:) = SDII_era_mean ;fi(1,:,:) = RX5d_era_mean ;fi(2,:,:) = R95pTOT_era_mean ; ;type = "pdf" ; choose x11 for pop-up window, or pdf, or eps ;file_save = "Plot/R95pTOT_era" + domain ;tiMainString = (/"W5k","W5k","W5k"/) ;txString = (/"(a) SDII","(b) RX5d","(c) R95pTOT"/) ;panel_x = 1 ;panel_y = 1 ;cnMinLevelValF = 0 ;cnMaxLevelValF = 250 ;cnLevelSpacingF =10 ;lbTitleString = (/"mm"/) ;plot = plot_cordex_panel_scale(cen_lat,cen_lon,type,file_save,fi,lat,lon,tiMainString,lbTitleString,panel_x,panel_y,cnMinLevelValF,cnMaxLevelValF,cnLevelSpacingF) ; ;delete(fi) ; ;fi = new((/1,x,y/),float) ; ;fi(0,:,:) = R99pTOT_era_mean ; ;type = "pdf" ; choose x11 for pop-up window, or pdf, or eps ;file_save = "Plot/R99pTOT_era" + domain ;tiMainString = (/"R99pTOT"/) ;panel_x = 1 ;panel_y = 1 ;cnMinLevelValF = 0 ;cnMaxLevelValF = 80 ;cnLevelSpacingF = 10 ;lbTitleString = (/"mm"/) ;plot = plot_cordex_panel_scale(cen_lat,cen_lon,type,file_save,fi,lat,lon,tiMainString,lbTitleString,panel_x,panel_y,cnMinLevelValF,cnMaxLevelValF,cnLevelSpacingF) ; ;delete(fi) ; ;fi = new((/1,x,y/),float) ; ;fi(0,:,:) = PRCPTOT_era_mean_mask ; ;type = "pdf" ; choose x11 for pop-up window, or pdf, or eps ;file_save = "Plot/PRCPTOT_era" + domain ;tiMainString = (/"PRCPTOT"/) ;panel_x = 1 ;panel_y = 1 ;cnMinLevelValF = 0 ;cnMaxLevelValF = 150 ;cnLevelSpacingF = 10 ;lbTitleString = (/"days"/) ;plot = plot_cordex_panel_scale(cen_lat,cen_lon,type,file_save,fi,lat,lon,tiMainString,lbTitleString,panel_x,panel_y,cnMinLevelValF,cnMaxLevelValF,cnLevelSpacingF) ; ;fi = new((/1,x,y/),float) ; ;fi(0,:,:) = PRCPTOT_bom_mean_mask ; ;type = "pdf" ; choose x11 for pop-up window, or pdf, or eps ;file_save = "Plot/PRCPTOT_bom" + domain ;tiMainString = (/"PRCPTOT"/) ;panel_x = 1 ;panel_y = 1 ;cnMinLevelValF = 0 ;cnMaxLevelValF = 150 ;cnLevelSpacingF = 10 ;lbTitleString = (/"days"/) ;plot = plot_cordex_panel_scale(cen_lat,cen_lon,type,file_save,fi,lat,lon,tiMainString,lbTitleString,panel_x,panel_y,cnMinLevelValF,cnMaxLevelValF,cnLevelSpacingF) ; ;fi = new((/1,x,y/),float) ; ;fi(0,:,:) = SPII_era_mean ; ;type = "pdf" ; choose x11 for pop-up window, or pdf, or eps ;file_save = "Plot/SPII_era" + domain ;tiMainString = (/"SPII"/) ;panel_x = 1 ;panel_y = 1 ;cnMinLevelValF = 0 ;cnMaxLevelValF = 12 ;cnLevelSpacingF = 1 ;lbTitleString = (/"mm/day"/) ;plot = plot_cordex_panel_scale(cen_lat,cen_lon,type,file_save,fi,lat,lon,tiMainString,lbTitleString,panel_x,panel_y,cnMinLevelValF,cnMaxLevelValF,cnLevelSpacingF) ; end