; 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 = "1997" days_in_year = 365 days_in_feb = 28 name = "BOM" diro = "Processed/" ;wrf inputs wrf_inp = (/"wrfhrly_d03_1993-01-01_00:00:00","wrfhrly_d02_1995-10-01_00:00:00"/) ;dimensions of domain x = (/691,691/) y = (/886,886/) filo_d03 = name + "_d03_post_regrid_" + year + ".nc" filo_d02 = name + "_d02_post_regrid" + year + ".nc" filo = (/filo_d03, filo_d02/) ;-------For each wrf domain ------------------ do i = 0,dimsizes(x)-1,1 ;------ Extract Monthly Values --------------- ; preallocate arrays for monthly data TXx = new((/12,x(i),y(i)/),float) TXn = new((/12,x(i),y(i)/),float) TNx = new((/12,x(i),y(i)/),float) TNn = new((/12,x(i),y(i)/),float) DTR = new((/12,x(i),y(i)/),float) TMAX = new((/12,x(i),y(i)/),float) TMIN = new((/12,x(i),y(i)/),float) FD = new((/12,x(i),y(i)/),float) SU = new((/12,x(i),y(i)/),float) TR = new((/12,x(i),y(i)/),float) PPT = new((/12,x(i),y(i)/),float) RX1d = new((/12,x(i),y(i)/),float) RX5d = new((/12,x(i),y(i)/),float) R10mm = new((/12,x(i),y(i)/),float) R20mm = new((/12,x(i),y(i)/),float) DLY_PPT = new((/days_in_year,x(i),y(i)/),float) DLY_TMAX = new((/days_in_year,x(i),y(i)/),float) DLY_TMIN = new((/days_in_year,x(i),y(i)/),float) ; open wrf data files for lat long data to regrid file_open_wrf = addfile(wrf_inp(i) + ".nc","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_TMAX = 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_TMIN = f_tmin[:]->data_tmin delete(files_tmin) f_ppt = addfiles(files_ppt + ".nc", "r") ListSetType (f_ppt, "join") dly_ppt = f_ppt[:]->data_ppt delete(files_ppt) delete(f_tmax) delete(f_tmin) delete(f_ppt) dsize_dtmax = dimsizes(dly_TMAX) print(dsize_dtmax) dsize_lat = dimsizes(bom_lat) print(dsize_lat) dsize_lon = dimsizes(bom_lon) print(dsize_lon) ;regrid onto wrf domain ;now compute the indices dependant on TMAX and TMIN TXx(counter,:,:) = dim_max_n(dly_TMAX,0) TXn(counter,:,:) = dim_min_n(dly_TMAX,0) TNx(counter,:,:) = dim_max_n(dly_TMIN,0) TNn(counter,:,:) = dim_min_n(dly_TMIN,0) DTR_m = dly_TMAX - dly_TMIN DTR(counter,:,:) = dim_avg_n(DTR_m,0) TMAX(counter,:,:) = dim_avg_n(dly_TMAX,0) TMIN(counter,:,:) = dim_avg_n(dly_TMIN,0) delete(DTR_m) FD(counter,:,:) = dim_num_n(dly_TMIN.lt.0,0) SU(counter,:,:) = dim_num_n(dly_TMAX.gt.25,0) TR(counter,:,:) = dim_num_n(dly_TMIN.gt.20,0) DLY_TMAX(start_ind:end_ind,:,:) = dly_TMAX DLY_TMIN(start_ind:end_ind,:,:) = dly_TMIN delete(dly_TMAX) delete(dly_TMIN) ;now compute the indices dependant on daily ppt PPT(counter,:,:) = dim_sum_n(dly_ppt,0) R10mm(counter,:,:) = dim_num_n(dly_ppt.gt.10,0) R20mm(counter,:,:) = dim_num_n(dly_ppt.gt.20,0) RX1d(counter,:,:) = dim_max_n(dly_ppt,0) ppt_rave = runave_n(dly_ppt,5,0,0) ppt_rave_max = dim_max_n(ppt_rave,0) RX5d(counter,:,:) = ppt_rave_max * 5 delete(ppt_rave) delete(ppt_rave_max) DLY_PPT(start_ind:end_ind,:,:) = dly_ppt delete(dly_ppt) start_ind = start_ind + days(counter) counter = counter + 1 end_ind = end_ind + days(counter) end do ; now that we have our years worth of data ; regrid ;regrid onto wrf domain ;first flip arrays size_bom_lat = dimsizes(bom_lat) TXx_flip = TXx(:,(size_bom_lat-1):0,:) TXn_flip = TXn(:,(size_bom_lat-1):0,:) TNx_flip = TNx(:,(size_bom_lat-1):0,:) TNn_flip = TNn(:,(size_bom_lat-1):0,:) DTR_flip = DTR(:,(size_bom_lat-1):0,:) TMIN_flip = TMIN(:,(size_bom_lat-1):0,:) TMAX_flip = TMAX(:,(size_bom_lat-1):0,:) PPT_flip = PPT(:,(size_bom_lat-1):0,:) FD_flip = FD(:,(size_bom_lat-1):0,:) SU_flip = SU(:,(size_bom_lat-1):0,:) TR_flip = TR(:,(size_bom_lat-1):0,:) RX1d_flip = RX1d(:,(size_bom_lat-1):0,:) RX5d_flip = RX5d(:,(size_bom_lat-1):0,:) R10mm_flip = R10mm(:,(size_bom_lat-1):0,:) R20mm_flip = R20mm(:,(size_bom_lat-1):0,:) DLY_PPT_flip = DLY_PPT(:,(size_bom_lat-1):0,:) DLY_TMAX_flip = DLY_TMAX(:,(size_bom_lat-1):0,:) DLY_TMIN_flip = DLY_TMIN(:,(size_bom_lat-1):0,:) delete(TXx) delete(TXn) delete(TNx) delete(TNn) delete(DTR) delete(TMAX) delete(TMIN) delete(FD) delete(SU) delete(TR) delete(PPT) delete(RX1d) delete(RX5d) delete(DLY_PPT) delete(DLY_TMAX) delete(DLY_TMIN) delete(R10mm) delete(R20mm) TXx = rgrid2rcm(bom_lat,bom_lon,TXx_flip,wrf_lat,wrf_lon,0) TXn = rgrid2rcm(bom_lat,bom_lon,TXn_flip,wrf_lat,wrf_lon,0) TNx = rgrid2rcm(bom_lat,bom_lon,TNx_flip,wrf_lat,wrf_lon,0) TNn = rgrid2rcm(bom_lat,bom_lon,TNn_flip,wrf_lat,wrf_lon,0) DTR = rgrid2rcm(bom_lat,bom_lon,DTR_flip,wrf_lat,wrf_lon,0) TMAX = rgrid2rcm(bom_lat,bom_lon,TMAX_flip,wrf_lat,wrf_lon,0) TMIN = rgrid2rcm(bom_lat,bom_lon,TMIN_flip,wrf_lat,wrf_lon,0) FD = rgrid2rcm(bom_lat,bom_lon,FD_flip,wrf_lat,wrf_lon,0) SU = rgrid2rcm(bom_lat,bom_lon,SU_flip,wrf_lat,wrf_lon,0) TR = rgrid2rcm(bom_lat,bom_lon,TR_flip,wrf_lat,wrf_lon,0) PPT = rgrid2rcm(bom_lat,bom_lon,PPT_flip,wrf_lat,wrf_lon,0) RX1d = rgrid2rcm(bom_lat,bom_lon,RX1d_flip,wrf_lat,wrf_lon,0) RX5d = rgrid2rcm(bom_lat,bom_lon,RX1d_flip,wrf_lat,wrf_lon,0) R20mm = rgrid2rcm(bom_lat,bom_lon,R20mm_flip,wrf_lat,wrf_lon,0) R10mm = rgrid2rcm(bom_lat,bom_lon,R10mm_flip,wrf_lat,wrf_lon,0) DLY_TMIN = rgrid2rcm(bom_lat,bom_lon,DLY_TMIN_flip,wrf_lat,wrf_lon,0) DLY_TMAX = rgrid2rcm(bom_lat,bom_lon,DLY_TMAX_flip,wrf_lat,wrf_lon,0) DLY_PPT = rgrid2rcm(bom_lat,bom_lon,DLY_PPT_flip,wrf_lat,wrf_lon,0) delete(TXx_flip) delete(TXn_flip) delete(TNx_flip) delete(TNn_flip) delete(DTR_flip) delete(TMAX_flip) delete(TMIN_flip) delete(FD_flip) delete(SU_flip) delete(TR_flip) delete(PPT_flip) delete(RX1d_flip) delete(RX5d_flip) delete(DLY_PPT_flip) delete(DLY_TMAX_flip) delete(DLY_TMIN_flip) delete(R10mm_flip) delete(R20mm_flip) ;---- write to a netcdf file for further processing --------- system("/bin/rm -f " + filo(i)) fout = addfile(diro+filo(i),"c") fout->TXx = (/TXx/) fout->TXn = (/TXn/) fout->TNx = (/TNx/) fout->TNn = (/TNn/) fout->DTR = (/DTR/) fout->TMAX_mean = (/TMAX/) fout->TMIN_mean = (/TMIN/) fout->FD = (/FD/) fout->SU = (/SU/) fout->TR = (/TR/) fout->PPT_mthly = (/PPT/) fout->RX1d = (/RX1d/) fout->RX5d = (/RX5d/) fout->DLY_PPT = (/DLY_PPT/) fout->DLY_TMAX = (/DLY_TMAX/) fout->DLY_TMIN = (/DLY_TMIN/) fout->R10mm = (/R10mm/) fout->R20mm = (/R20mm/) delete(TXx) delete(TXn) delete(TNx) delete(TNn) delete(DTR) delete(TMAX) delete(TMIN) delete(FD) delete(SU) delete(TR) delete(PPT) delete(RX1d) delete(RX5d) delete(DLY_PPT) delete(DLY_TMAX) delete(DLY_TMIN) delete(R10mm) delete(R20mm) delete(file_open_wrf) delete(wrf_lat) delete(wrf_lon) end do end