; Julia 30/08/2013: script to process and regrid BOM-AWAP ; files for precipitation and temperature metrics and ; indices. This script is designed to be run in the ; directory containing the dirs for BOM PPT TMAX and TMIN ; files 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" begin ;****** INPUTS ************ year_start = 1970 year_end = 2010 name = "BOM_lambert" diro = "Processed/" ;wrf inputs wrf_inp = "wrfout_d02.nc" ;dimensions of domain x = 140 y = 175 ;************************* ;loop through years years = ispan(year_start,year_end,1) years@calendar = "standard" ;years = (/1970,1971,1972,1973,1974,1975,1976,1977,1978,1979,1980,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,2010/) Days_in_year = new((/dimsizes(years)/),float) Days_in_feb = new((/dimsizes(years)/),float) isleap = isleapyear(years) Days_in_year = where(isleap.eq.True,366,365) Days_in_year = where(isleap.eq.True,29,28) ;Days_in_feb = (/28,28,29,28,28,28,29,28,28,28,29,28,28,28,29,28,28,28,29,28,28,28,29,28,28,28,29,28,28,28,29,28,28,28,29,28,28,28,29,28,28/) ;Days_in_year = (/365,365,366,365,365,365,366,365,365,365,366,365,365,365,366,365,365,365,366,365,365,365,366,365,365,365,366,365,365,365,366,365,365,365,366,365,365,365,366,365,365/) do i = 0,dimsizes(years)-1,1 year = years(i) days_in_year = Days_in_year(i) days_in_feb = Days_in_feb(i) filo = name + "_d02_" + year + ".nc" ;------ Extract Monthly Values --------------- ; preallocate arrays for monthly data TMAX = new((/12,x,y/),float) TMIN = new((/12,x,y/),float) PPT = new((/12,x,y/),float) ; open wrf data files for lat long data to regrid file_open_wrf = addfile(wrf_inp,"r") wrf_lat = file_open_wrf->XLAT(0,:,:) wrf_lon = file_open_wrf->XLONG(0,:,:) ; now get data a month at a time mm = (/1,2,3,4,5,6,7,8,9,10,11,12/) months = (/"01","02","03","04","05","06","07","08","09","10","11","12"/) days = (/31,days_in_feb,31,30,31,30,31,31,30,31,30,31,31/) print(months) print(days) counter = 0; start_ind = 0 end_ind = days(0)-1 do j = 0,dimsizes(months)-1,1 print(end_ind) print(start_ind) month = months(j) files_ppt = systemfunc("ls ../ppt/" + year + month + "*" + year + month + "*") files_tmax = systemfunc("ls ../tmax/" + year + month + "*" + year + month + "*") files_tmin = systemfunc("ls ../tmin/" + year + month + "*" + year + month + "*") f_tmax = addfiles(files_tmax + ".nc", "r") ListSetType (f_tmax, "join") dly_MAX = f_tmax[:]->data_tmax bom_lat = f_tmax[0]->lat_tmax bom_lon = f_tmax[0]->lon_tmax delete(files_tmax) f_tmin = addfiles(files_tmin + ".nc", "r") ListSetType (f_tmin, "join") dly_MIN = f_tmin[:]->data_tmin delete(files_tmin) f_ppt = addfiles(files_ppt + ".nc", "r") ListSetType (f_ppt, "join") dly_prec = f_ppt[:]->data_ppt delete(files_ppt) delete(f_tmax) delete(f_tmin) delete(f_ppt) dsize_dtmax = dimsizes(dly_MAX) print(dsize_dtmax) dsize_lat = dimsizes(bom_lat) print(dsize_lat) dsize_lon = dimsizes(bom_lon) print(dsize_lon) ;regrid onto wrf domain ;first flip bom array size_bom_lat = dimsizes(bom_lat) dly_TMAX_flip = dly_MAX(:,(size_bom_lat-1):0,:) dly_TMIN_flip = dly_MIN(:,(size_bom_lat-1):0,:) dly_ppt_flip = dly_prec(:,(size_bom_lat-1):0,:) delete(dly_MAX) delete(dly_MIN) delete(dly_prec) dsize_tmax_flip = dimsizes(dly_TMAX_flip) print(dsize_tmax_flip) dly_TMAX = rgrid2rcm(bom_lat,bom_lon,dly_TMAX_flip,wrf_lat,wrf_lon,0) dly_TMIN = rgrid2rcm(bom_lat,bom_lon,dly_TMIN_flip,wrf_lat,wrf_lon,0) dly_ppt = rgrid2rcm(bom_lat,bom_lon,dly_ppt_flip,wrf_lat,wrf_lon,0) delete(dly_TMAX_flip) delete(dly_TMIN_flip) delete(dly_ppt_flip) dsizes_dly_TMAX = dimsizes(dly_TMAX) print(dsizes_dly_TMAX) TMAX(counter,:,:) = dim_avg_n(dly_TMAX,0) TMIN(counter,:,:) = dim_avg_n(dly_TMIN,0) PPT(counter,:,:) = dim_sum_n(dly_ppt,0) delete(dly_TMAX) delete(dly_TMIN) delete(dly_ppt) counter = counter + 1 end do ; now that we have our years worth of data ;---- write to a netcdf file for further processing --------- system("/bin/rm -f " + filo) fout = addfile(diro+filo,"c") fout->TMAX_mean = (/TMAX/) fout->TMIN_mean = (/TMIN/) fout->PPT_mthly = (/PPT/) delete(TMAX) delete(TMIN) delete(PPT) delete(file_open_wrf) delete(wrf_lat) delete(wrf_lon) delete(counter) end do end