; Julia 26/08/2013: script to process wrfhrly files for ; precipitation and temperature metrics and indices. This ; script is designed to be run in the directory where the ; wrfhrly files reside. Script also adjusts wrf UTC output onto ; the BOM daily timescale. (BOM data is UTC +9 and days run ; from 9am to 9am) 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 = 1981 year_1 = year + 1 days_in_year = 365 days_in_feb = 28 name = "ERA" domain = "d03" file_type = "wrfhrly" diro = "Processed/" ;dimensions of domain, less relaxation points (10 on each side) x = 100 y = 100 filo = name + "_" + domain + "_" + year + "_tzadjust.nc" ;------ Extract Monthly Values --------------- ; preallocate arrays for monthly data TXx = new((/12,x,y/),float) TXn = new((/12,x,y/),float) TNx = new((/12,x,y/),float) TNn = new((/12,x,y/),float) DTR = new((/12,x,y/),float) TMAX = new((/12,x,y/),float) TMIN = new((/12,x,y/),float) FD = new((/12,x,y/),float) SU = new((/12,x,y/),float) TR = new((/12,x,y/),float) PPT = new((/12,x,y/),float) RX1d = new((/12,x,y/),float) RX5d = new((/12,x,y/),float) R10mm = new((/12,x,y/),float) R20mm = new((/12,x,y/),float) DLY_PPT = new((/days_in_year,x,y/),float) DLY_TMAX = new((/days_in_year,x,y/),float) DLY_TMIN = new((/days_in_year,x,y/),float) ; open a wrfout file to extract the land mask fi = addfile("wrfout_" + domain + ".nc", "r") landmask = wrf_user_getvar(fi,"LANDMASK",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","01"/) 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 i = 0,dimsizes(mm)-1,1 print(end_ind) print(start_ind) month = months(i) month_1 = months(counter+1) files_1 = systemfunc("ls " + file_type + "_" + domain + "_" + year + "-" + month + "*") files_2 = systemfunc("ls " + file_type + "_" + domain + "_" + year + "-" + month_1 + "-01_00:00:00") files_3 = systemfunc("ls " + file_type + "_" + domain + "_" + year_1 + "-" + month_1 + "-01_00:00:00") if ((month).eq.12) then files = array_append_record(files_1, files_3,0) else files = array_append_record(files_1, files_2,0) end if print(files) f = addfiles(files + ".nc", "r") ListSetType (f, "cat") T = f[:]->T2 T_c = T - 273.15 T_mask = mask(T_c,landmask,1) dsizes_array = dimsizes(T_mask) ;Remove 10 relaxation points on all sides T_norelax = T_mask(:,10:(dsizes_array(1)-11),10:(dsizes_array(2)-11)) delete(files_1) delete(files_2) delete(files) delete(T) delete(T_c) delete(T_mask) ;preallocate array for daily min and max temps dly_TMAX = new((/days(counter),x,y/),"float") dly_TMIN = new((/days(counter),x,y/),"float") ;loop through every 24 hours and extract max and min temps size_temp_array_T2 = dimsizes(T_norelax) print(size_temp_array_T2) hour_index_start_T2 = ispan(1,(size_temp_array_T2(0)-24),24) hour_index_end_T2 = ispan(24,(size_temp_array_T2(0)-24),24) print(hour_index_start_T2) print(hour_index_end_T2) do j = 0,dimsizes(hour_index_start_T2)-1,1 int_24_temps = T_norelax(hour_index_start_T2(j):hour_index_end_T2(j),:,:) wrf_max_temps = dim_max_n(int_24_temps,0) wrf_min_temps = dim_min_n(int_24_temps,0) dly_TMAX(j,:,:) = wrf_max_temps dly_TMIN(j,:,:) = wrf_min_temps delete(int_24_temps) delete(wrf_max_temps) delete(wrf_min_temps) end do delete(size_temp_array_T2) delete(hour_index_start_T2) delete(hour_index_end_T2) ;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) delete(T_norelax) ;we are now extracting rainfall metrics NC = f[:]->RAINNC RAINC = f[:]->RAINC ppt = NC + RAINC ;now mask out land and remove 10 relaxation points on all sides ppt_mask = mask(ppt,landmask,1) ppt_norelax = ppt_mask(:,10:(dsizes_array(1)-11),10:(dsizes_array(2)-11)) delete(NC) delete(RAINC) delete(ppt) delete(ppt_mask) ;preallocate array for daily ppt dly_ppt = new((/days(counter),x,y/),"float") ;loop through 24 hours and extract daily rainfall size_temp_array_ppt = dimsizes(ppt_norelax) hour_index_start_ppt = ispan(1,(size_temp_array_ppt(0)-24),24) hour_index_end_ppt = ispan(24,(size_temp_array_ppt(0)-24),24) do k = 0,dimsizes(hour_index_start_ppt)-1,1 int_24_ppt_start = ppt_norelax(hour_index_start_ppt(k),:,:) int_24_ppt_end = ppt_norelax(hour_index_end_ppt(k),:,:) sum_ppt = int_24_ppt_end - int_24_ppt_start dly_ppt(k,:,:) = sum_ppt delete(int_24_ppt_start) delete(int_24_ppt_end) delete(sum_ppt) end do delete(size_temp_array_ppt) delete(hour_index_start_ppt) delete(hour_index_end_ppt) delete(ppt_norelax) ;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 ;---- write to a netcdf file for further processing --------- system("/bin/rm -f " + filo) fout = addfile(diro+filo,"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) end