; Julia 4/11/2013: this script loads in the yearly data files and extracts the daily data from these files and then assigns these values in a probability density function. The pdf of the model and observatio are then compared to produce a perkins skill score at each grid point. These are then plotted on a contour. 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" begin ;---------------INPUTS------------------------------------- name = "ERA" file_suffix = "tzadjust" domain = "d02" filo = name + "_" + domain + "_pkss.nc" diro = "Processed/" diro_bom = "/group/y98/julia/BOM/Processed/" files_input_era = systemfunc("ls " + diro + name + "_" + domain + "_*_tzadjust.nc") files_input_bom = systemfunc("ls " + diro_bom + "BOM_" + domain + "_*.nc") x = 125 y = 160 ;------load yearly files --------------- f_era = addfiles(files_input_era + ".nc", "r") ListSetType (f_era, "cat") f_bom = addfiles(files_input_bom + ".nc", "r") ListSetType (f_bom, "cat") ; -------- now extract the vairables we need TMAX_era = f_era[:]->DLY_TMAX TMIN_era = f_era[:]->DLY_TMIN PPT_era = f_era[:]->DLY_PPT TMAX_bom = f_bom[:]->DLY_TMAX TMIN_bom = f_bom[:]->DLY_TMIN PPT_bom = f_bom[:]->DLY_PPT printVarSummary(TMAX_bom) printVarSummary(TMAX_era) 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) ; we need to remove the relaxation zone from the BOM data too and we have one extra day of ERA data that I need to remove dsizes_era = dimsizes(TMAX_era) TMAX_era_adj = TMAX_era(1:(dsizes_era(0)-1),:,:) TMIN_era_adj = TMIN_era(1:(dsizes_era(0)-1),:,:) PPT_era_adj = PPT_era(1:(dsizes_era(0)-1),:,:) 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)) PPT_bom_nr = PPT_bom(:,10:(dsizes_array(0)-11),10:(dsizes_array(1)-11)) tmax_era_mask = mask(TMAX_era_adj,landmask_nr,1) tmin_era_mask = mask(TMIN_era_adj,landmask_nr,1) ppt_era_mask = mask(PPT_era_adj,landmask_nr,1) tmax_bom_mask = mask(TMAX_bom_nr,landmask_nr,1) tmin_bom_mask = mask(TMIN_bom_nr,landmask_nr,1) ppt_bom_mask = mask(PPT_bom_nr,landmask_nr,1) delete(TMAX_era) delete(TMIN_era) delete(PPT_era) delete(TMAX_bom) delete(TMIN_bom) delete(PPT_bom) delete(PPT_bom_nr) delete(TMAX_bom_nr) delete(TMIN_bom_nr) delete(PPT_era_adj) delete(TMAX_era_adj) delete(TMIN_era_adj) ppt_era_mask_adj = where(ppt_era_mask.lt.(0.2),(9.96921e+36),ppt_era_mask) delete(ppt_era_mask) ppt_bom_mask_adj = where(ppt_bom_mask.lt.(0.2),(9.96921e+36),ppt_bom_mask) delete(ppt_bom_mask) ppt_era_mask_adj@_FillValue = 9.96921e+36 ppt_bom_mask_adj@_FillValue = 9.96921e+36 ppt_pkss = new((/x,y/),float) tmax_pkss = new((/x,y/),float) tmin_pkss = new((/x,y/),float) counteri = 0 do i = 0,x-1,1 counterj = 0 do j = 0,y-1,1 ;***** tmax ****** opt = True opt@bin_min = 9. opt@bin_max = 45. tmax_era_pdfx = pdfx(tmax_era_mask(:,counteri,counterj),36,opt) tmax_bom_pdfx = pdfx(tmax_bom_mask(:,counteri,counterj),36,opt) tmax_era_pdfx_float = doubletofloat(tmax_era_pdfx/100) tmax_bom_pdfx_float = doubletofloat(tmax_bom_pdfx/100) tmax_perkins = new((/36/),float) do l = 0,35,1 tmax_perkins_temp = new((/2/),float) tmax_perkins_temp(0) = tmax_era_pdfx_float(l) tmax_perkins_temp(1) = tmax_bom_pdfx_float(l) tmax_perkins(l) = min(tmax_perkins_temp) delete(tmax_perkins_temp) end do sum_tmax_perkins = sum(tmax_perkins) tmax_pkss(counteri,counterj) = sum_tmax_perkins ;delete(opt@bin_min) ;delete(opt@bin_max) ;delete(tmax_era_pdfx) ;delete(tmax_bom_pdfx) ;delete(tmax_era_pdfx_float) ;delete(tmax_bom_pdfx_float) ;delete(tmax_perkins) ;delete(sum_tmax_perkins) ;***** tmin ****** opt@bin_min = -1 opt@bin_max = 30 tmin_era_pdfx = pdfx(tmin_era_mask(:,counteri,counterj),31,opt) tmin_bom_pdfx = pdfx(tmin_bom_mask(:,counteri,counterj),31,opt) tmin_era_pdfx_float = doubletofloat(tmin_era_pdfx/100) tmin_bom_pdfx_float = doubletofloat(tmin_bom_pdfx/100) tmin_perkins = new((/31/),float) do m = 0,30,1 tmin_perkins_temp = new((/2/),float) tmin_perkins_temp(0) = tmin_era_pdfx_float(m) tmin_perkins_temp(1) = tmin_bom_pdfx_float(m) tmin_perkins(m) = min(tmin_perkins_temp) ;delete(tmin_perkins_temp) end do sum_tmin_perkins = sum(tmin_perkins) tmin_pkss(counteri,counterj) = sum_tmin_perkins ;delete(opt@bin_min) ;delete(opt@bin_max) ;delete(tmin_era_pdfx) ;delete(tmin_bom_pdfx) ;delete(tmin_era_pdfx_float) ;delete(tmin_bom_pdfx_float) ;delete(tmin_perkins) ;delete(sum_tmin_perkins) ; ******** ppt ******** opt@bin_min = 0 opt@bin_max = 20 ppt_era_pdfx = pdfx(ppt_era_mask_adj(:,counteri,counterj),20,opt) ppt_bom_pdfx = pdfx(ppt_bom_mask_adj(:,counteri,counterj),20,opt) ppt_era_pdfx_float = doubletofloat(ppt_era_pdfx/100) ppt_bom_pdfx_float = doubletofloat(ppt_bom_pdfx/100) ppt_perkins = new((/20/),float) do k = 0,19,1 ppt_perkins_temp = new((/2/),float) ppt_perkins_temp(0) = ppt_era_pdfx_float(k) ppt_perkins_temp(1) = ppt_bom_pdfx_float(k) ppt_perkins(k) = min(ppt_perkins_temp) ;delete(ppt_perkins_temp) end do sum_ppt_perkins = sum(ppt_perkins) ;print(sum_ppt_perkins) ppt_pkss(counteri,counterj) = sum_ppt_perkins ;delete(opt@bin_min) ;delete(opt@bin_max) ;delete(ppt_era_pdfx) ;delete(ppt_bom_pdfx) ;delete(ppt_era_pdfx_float) ;delete(ppt_bom_pdfx_float) ;delete(ppt_perkins) ;delete(sum_ppt_perkins) counterj = counterj + 1 end do counteri = counteri + 1 end do ; we need to limit the d02 grid to the extent of the d03. We do this using predetermined ij locations(see plot_contours_func_tave_d02.ncl). ppt_pkss_limit = ppt_pkss(45:94,71:120) tmax_pkss_limit = tmax_pkss(45:94,71:120) tmin_pkss_limit = tmin_pkss(45:94,71:120) lat_limit = lat(45:94,71:120) lon_limit = lon(45:94,71:120) dimsizes_limit = dimsizes(tmax_pkss_limit) print("maximum ppt_pkss is") print(max(ppt_pkss_limit)) print("maximum tmax_pkss is") print(max(tmax_pkss_limit)) print("maximum tmin_pkss is") print(max(tmin_pkss_limit)) print("minimum ppt_pkss is") print(min(ppt_pkss_limit)) print("minimum tmax_pkss is") print(min(tmax_pkss_limit)) print("minimum tmin_pkss is") print(min(tmin_pkss_limit)) print("mean ppt_pkss is") print(avg(ppt_pkss_limit)) print("mean tmax_pkss is") print(avg(tmax_pkss_limit)) print("mean tmin_pkss is") print(avg(tmin_pkss_limit)) ;---- write to a netcdf file for further processing --------- system("/bin/rm -f " + filo) fout = addfile(diro+filo,"c") fout->tmin_pkss_d02 = (/tmin_pkss_limit/) fout->tmax_pkss_d02 = (/tmax_pkss_limit/) fout->ppt_pkss_d02 = (/ppt_pkss_limit/) ;*******Draw Panel Plots************ printVarSummary(ppt_pkss) printVarSummary(ppt_pkss_limit) printVarSummary(lat) printVarSummary(lat_limit) fi = new((/1,dimsizes_limit(0),dimsizes_limit(1)/),float) fi(0,:,:) = ppt_pkss_limit(:,:) ;fi(1,:,:) = tmax_pkss_limit(:,:) ;fi(2,:,:) = tmin_pkss_limit(:,:) type = "pdf" ; choose x11 for pop-up window, or pdf, or eps file_save = "Plot/Perkins_score_ppt_" + domain tiMainString = (/""/) panel_x = 1 panel_y = 1 txString = "(a) W10k" lbTitleString = "" cnMinLevelValF = 0.725 cnMaxLevelValF = 0.975 cnLevelSpacingF = 0.025 plot = plot_cordex_panel_scale(cen_lat,cen_lon,type,file_save,fi,lat_limit,lon_limit,txString,tiMainString,lbTitleString,panel_x,panel_y,cnMinLevelValF,cnMaxLevelValF,cnLevelSpacingF) delete(fi) delete(tiMainString) fi = new((/2,dimsizes_limit(0),dimsizes_limit(1)/),float) fi(0,:,:) = tmax_pkss_limit(:,:) fi(1,:,:) = tmin_pkss_limit(:,:) ;fi(2,:,:) = tmin_pkss_limit(:,:) type = "pdf" ; choose x11 for pop-up window, or pdf, or eps file_save = "Plot/Perkins_score_temp_" + domain tiMainString = (/"MAX","MIN"/) panel_x = 1 panel_y = 2 txString = "(a) W10k" lbTitleString = "" cnMinLevelValF = 0.725 cnMaxLevelValF = 0.975 cnLevelSpacingF = 0.025 plot = plot_cordex_panel_scale(cen_lat,cen_lon,type,file_save,fi,lat_limit,lon_limit,txString,tiMainString,lbTitleString,panel_x,panel_y,cnMinLevelValF,cnMaxLevelValF,cnLevelSpacingF) end