; Julia 11/11/2013: this script loads in the daily data from the yearly data files. We are calculating the rx5day by calculating the 5 day running mean. then we work out the max 5 days running mean and multiply this by 5, giving us the maximum 5 day precipitation rate. 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------------------------------------- FD_era = new((/29,12,100,100/),float) FD_bom = new((/29,12,100,100/),float) SU_era = new((/29,12,100,100/),float) SU_bom = new((/29,12,100,100/),float) do i = 0,28,1 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/) months = (/1,2,3,4,5,6,7,8,9,10,11,12/) name = "ERA" file_suffix = "tzadjust" domain = "d03" filo = name + "_" + domain + "_temp_indices_gs.nc" diro = "Processed/" diro_bom = "/group/y98/julia/BOM/Processed/" files_input_era = systemfunc("ls " + diro + name + "_" + domain + "_" + year(i) + "_tzadjust.nc") files_input_bom = systemfunc("ls " + diro_bom + "BOM_" + domain + "_" +year(i) + ".nc") x = 100 y = 100 ;------load files, a year at a time and work out Rx5d, CWD, CDD --------------- 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 tmax_era = f_era[:]->DLY_TMAX tmax_bom = f_bom[:]->DLY_TMAX tmin_era = f_era[:]->DLY_TMIN tmin_bom = f_bom[:]->DLY_TMIN delete(f_bom) delete(f_era) ;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) 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) print(year(i)) ; we need to remove the relaxation zone from the BOM data too tmax_bom_nr = tmax_bom(:,10:(dsizes_array(0)-11),10:(dsizes_array(1)-11)) tmin_bom_nr = tmin_bom(:,10:(dsizes_array(0)-11),10:(dsizes_array(1)-11)) tmax_era_mask = mask(tmax_era,landmask_nr,1) tmin_era_mask = mask(tmin_era,landmask_nr,1) tmax_bom_mask = mask(tmax_bom_nr,landmask_nr,1) tmin_bom_mask = mask(tmin_bom_nr,landmask_nr,1) delete(tmax_era) delete(tmax_bom) delete(tmax_bom_nr) delete(tmin_era) delete(tmin_bom) delete(tmin_bom_nr) ; now we need to put this data back into months dsizes_tmax = dimsizes(tmax_bom_mask) if (dsizes_tmax(0).eq.365) then month_start = (/0,31,59,90,120,151,181,212,243,273,304,334/) month_end = (/30,58,89,119,150,180,211,242,272,303,333,364/) else month_start = (/0,31,60,91,121,151,182,213,244,274,305,335/) month_end = (/30,59,90,120,151,181,212,243,273,304,334,365/) end if FD_bom_monthly = new((/12,x,y/),float) FD_era_monthly = new((/12,x,y/),float) SU_bom_monthly = new((/12,x,y/),float) SU_era_monthly = new((/12,x,y/),float) do j = 0,dimsizes(month_start)-1,1 start_ind = month_start(j) end_ind = month_end(j) print(j) tmax_bom_month = tmax_bom_mask(start_ind:end_ind,:,:) tmax_era_month = tmax_era_mask(start_ind:end_ind,:,:) tmin_bom_month = tmin_bom_mask(start_ind:end_ind,:,:) tmin_era_month = tmin_era_mask(start_ind:end_ind,:,:) FD_bom_monthly(j,:,:) = dim_num_n(tmin_bom_month.lt.2,0) FD_era_monthly(j,:,:) = dim_num_n(tmin_era_month.lt.2,0) SU_bom_monthly(j,:,:) = dim_num_n(tmax_bom_month.ge.34,0) SU_era_monthly(j,:,:) = dim_num_n(tmax_era_month.ge.34,0) delete(tmax_bom_month) delete(tmax_era_month) delete(tmin_bom_month) delete(tmin_era_month) end do delete(tmin_bom_mask) delete(tmax_bom_mask) delete(tmin_era_mask) delete(tmax_era_mask) FD_era(i,:,:,:) = FD_era_monthly FD_bom(i,:,:,:) = FD_bom_monthly SU_era(i,:,:,:) = SU_era_monthly SU_bom(i,:,:,:) = SU_bom_monthly delete(FD_era_monthly) delete(FD_bom_monthly) delete(SU_bom_monthly) delete(SU_era_monthly) end do FD_era_mean = dim_avg_n(FD_era,0) FD_bom_mean = dim_avg_n(FD_bom,0) FD_bias = (FD_era_mean - FD_bom_mean) SU_era_mean = dim_avg_n(SU_era,0) SU_bom_mean = dim_avg_n(SU_bom,0) SU_bias = (SU_era_mean - SU_bom_mean) FD_era_mean_mask = mask(FD_era_mean,landmask_nr,1) FD_bom_mean_mask = mask(FD_bom_mean,landmask_nr,1) SU_era_mean_mask = mask(SU_era_mean,landmask_nr,1) SU_bom_mean_mask = mask(SU_bom_mean,landmask_nr,1) SU_bom_mean_mask_max = dim_max_n(SU_bom_mean,(/1,2/)) SU_era_mean_mask_max = dim_max_n(SU_era_mean,(/1,2/)) print(SU_bom_mean_mask_max) print(SU_era_mean_mask_max) fi = new((/6,x,y/),float) fi(0,:,:) = FD_bom_mean_mask(4,:,:) fi(1,:,:) = FD_bom_mean_mask(5,:,:) fi(2,:,:) = FD_bom_mean_mask(6,:,:) fi(3,:,:) = FD_bom_mean_mask(7,:,:) fi(4,:,:) = FD_bom_mean_mask(8,:,:) fi(5,:,:) = FD_bom_mean_mask(9,:,:) type = "pdf" ; choose x11 for pop-up window, or pdf, or eps file_save = "Plot/FD_bom_growing_season_" + domain tiMainString = (/"05","06","07","08","09","10"/) txString = "(c) Observations" panel_x = 1 panel_y = 6 cnMinLevelValF = 1 cnMaxLevelValF = 10 cnLevelSpacingF = 1 lbTitleString = (/"(Days)"/) plot = plot_cordex_panel_scale(cen_lat,cen_lon,type,file_save,fi,lat,lon,txString,tiMainString,lbTitleString,panel_x,panel_y,cnMinLevelValF,cnMaxLevelValF,cnLevelSpacingF) fi = new((/6,x,y/),float) fi(0,:,:) = FD_era_mean_mask(4,:,:) fi(1,:,:) = FD_era_mean_mask(5,:,:) fi(2,:,:) = FD_era_mean_mask(6,:,:) fi(3,:,:) = FD_era_mean_mask(7,:,:) fi(4,:,:) = FD_era_mean_mask(8,:,:) fi(5,:,:) = FD_era_mean_mask(9,:,:) type = "pdf" ; choose x11 for pop-up window, or pdf, or eps file_save = "Plot/FD_era_growing_season_" + domain tiMainString = (/"05","06","07","08","09","10"/) txString = "(b) W5k" panel_x = 1 panel_y = 6 cnMinLevelValF = 1 cnMaxLevelValF = 10 cnLevelSpacingF = 1 lbTitleString = (/"Days"/) plot = plot_cordex_panel_scale(cen_lat,cen_lon,type,file_save,fi,lat,lon,txString,tiMainString,lbTitleString,panel_x,panel_y,cnMinLevelValF,cnMaxLevelValF,cnLevelSpacingF) delete(fi) fi = new((/7,x,y/),float) fi(0,:,:) = SU_bom_mean_mask(0,:,:) fi(1,:,:) = SU_bom_mean_mask(1,:,:) fi(2,:,:) = SU_bom_mean_mask(2,:,:) fi(3,:,:) = SU_bom_mean_mask(3,:,:) fi(4,:,:) = SU_bom_mean_mask(9,:,:) fi(5,:,:) = SU_bom_mean_mask(10,:,:) fi(6,:,:) = SU_bom_mean_mask(11,:,:) delete(tiMainString) delete(txString) delete(plot) delete(cnLevelSpacingF) type = "pdf" ; choose x11 for pop-up window, or pdf, or eps file_save = "Plot/SU_34_bom_growing_season_" + domain tiMainString = (/"01","02","03","04","10","11","12"/) txString = "(c) Observations" panel_x = 2 panel_y = 4 cnMinLevelValF = 1 cnMaxLevelValF = 19 cnLevelSpacingF = 1 lbTitleString = (/"(Days)"/) plot = plot_cordex_panel_scale(cen_lat,cen_lon,type,file_save,fi,lat,lon,txString,tiMainString,lbTitleString,panel_x,panel_y,cnMinLevelValF,cnMaxLevelValF,cnLevelSpacingF) fi = new((/7,x,y/),float) fi(0,:,:) = SU_era_mean_mask(0,:,:) fi(1,:,:) = SU_era_mean_mask(1,:,:) fi(2,:,:) = SU_era_mean_mask(2,:,:) fi(3,:,:) = SU_era_mean_mask(3,:,:) fi(4,:,:) = SU_era_mean_mask(9,:,:) fi(5,:,:) = SU_era_mean_mask(10,:,:) fi(6,:,:) = SU_era_mean_mask(11,:,:) delete(cnLevelSpacingF) type = "pdf" ; choose x11 for pop-up window, or pdf, or eps file_save = "Plot/SU_34_era_growing_season_" + domain tiMainString = (/"01","02","03","04","10","11","12"/) txString = "(b) W5k" panel_x = 2 panel_y = 4 cnMinLevelValF = 1 cnMaxLevelValF = 19 cnLevelSpacingF = 1 lbTitleString = (/"Days"/) plot = plot_cordex_panel_scale(cen_lat,cen_lon,type,file_save,fi,lat,lon,txString,tiMainString,lbTitleString,panel_x,panel_y,cnMinLevelValF,cnMaxLevelValF,cnLevelSpacingF) delete(fi) end