; --------------------------------------------------------------------------- ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Program Name: combine.ncl, ; ; ; ; National Water Model (NWM) WRF-hydro forcing engine is developed by ; ; National Center for Atmospheric Research (NCAR), under the sponsorship ; ; of National Water Center (NWC). ; ; ; ; Team Members: ; ; NCAR Staff: Linlin Pan, Wei Yu, and David Gochis ; ; NWC/OWP Staff: Brian Cosgrove, Zhengtao Cui, Cham Pham, and James Taft ; ; ; ; This is a ncl program to combine file. ; ; ; ; Input: hrrr file, gfs file, output file ; ; ; ; Output: combined file ; ; ; ; For non-fatal errors output is witten to $DATA/logs ; ; ; ; Author(s)/Contact(s): Linlin Pan, lpan@ucar.edu ; ; Origination Sept., 2015 ; ; ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; load "$USHnwm/scale/set_scale_factor.ncl" ; filename1="$USHnwm/scale/set_scale_factor.ncl" ; if(isfilepresent(filename1)) then ; loadscript(filename1) ; else ; print("Error: scale factor file doesn't exist. Can't continue.") ; exit ; end if begin hrrrFile = getenv("hrrrFile") ; HRRR data gfsFile = getenv("gfsFile") ; GFS data outFile = getenv("outFile") ; outFile data init_time = getenv("init_time") ; outFile data valid_time = getenv("valid_time") ; outFile data print(valid_time) indexFlag = getenv("INDEX") dtime_map = (/4,2,2,2/) init_dtime = str_split_by_length(init_time, dtime_map) valid_dtime = str_split_by_length(valid_time, dtime_map) if(indexFlag .eq. "true") then indexf = addfile("index.nc","r") index = (/indexf->index/) end if f1 = addfile(hrrrFile,"r") f2 = addfile(gfsFile,"r") print("combine f1 :" + hrrrFile) print("combine f2 :" + gfsFile) system("rm -f "+outFile) fout = addfile(outFile,"c") globalAtt = True globalAtt@NWM_version_number=nwm_version globalAtt@model_output_type="forcing" globalAtt@model_configuration="short_range" globalAtt@model_total_valid_times=18 globalAtt@model_output_valid_time = valid_dtime(0)+"-"+valid_dtime(1)+"-"+valid_dtime(2)+"_"+valid_dtime(3)+":00:00" globalAtt@model_initialization_time = init_dtime(0)+"-"+init_dtime(1)+"-"+init_dtime(2)+"_"+init_dtime(3)+":00:00" units="minutes since 1970-01-01 00:00:00 UTC" year=tointeger(valid_dtime(0)) month=tointeger(valid_dtime(1)) day=tointeger(valid_dtime(2)) hour=tointeger(valid_dtime(3)) time = tointeger(cd_inv_calendar(year,month,day,hour,0,0,units, 0)) yeari=tointeger(init_dtime(0)) monthi=tointeger(init_dtime(1)) dayi=tointeger(init_dtime(2)) houri=tointeger(init_dtime(3)) houre=houri+21 time_min = tointeger(cd_inv_calendar(yeari,monthi,dayi,houri,0,0,units, 0)) time_max = tointeger(cd_inv_calendar(yeari,monthi,dayi,houre,0,0,units, 0)) fileattdef( fout, globalAtt ) reference_time = doubletoint(cd_inv_calendar(yeari,monthi,dayi,houri,0,0,units, 0)) reference_time!0="reference_time" reference_time@long_name="model initialization time" reference_time@standard_name="forecast_reference_time" reference_time@units=units dd = dimsizes(f1->T2D) ndd=dimsizes(dd) print("dd="+dd+" ndd="+ndd) dimlat = dd(ndd-2) dimlon = dd(ndd-1) ; dimNames = (/"time", "y", "x","reference_time","nv"/) ; dimSizes = (/ -1 , dimlat, dimlon , 1, 2/) ; dimUnlim = (/ True , False, False,False,False/) dimNames = (/"time", "y", "x","reference_time"/) dimSizes = (/ -1 , dimlat, dimlon , 1/) dimUnlim = (/ True , False, False,False/) filedimdef(fout,dimNames,dimSizes,dimUnlim) ; time=new((/1/),integer) ; time=stringtointeger(time) time!0="time" time@long_name="valid output time" ; time@bounds="time_bounds" time@standard_name="time" time@valid_max=time_max time@valid_min=time_min time@units = units filevardef(fout, "time" ,typeof(time) ,getvardims(time)) filevarattdef(fout,"time" ,time) fout->time=(/time/) ; filevardef(fout, "reference_time" ,typeof(reference_time) ,getvardims(reference_time)) filevarattdef(fout,"reference_time" ,reference_time) fout->reference_time=(/reference_time/) ; time_bounds=new((/1,2/),integer) time_bounds!0="time" time_bounds!1="nv" time_bounds@units=units time_bounds(0,0)=(/reference_time(0)/) time_bounds(0,1)=doubletoint(time(0)) delete(time_bounds@_FillValue) ; filevardef(fout, "time_bounds" ,typeof(time_bounds) ,getvardims(time_bounds)) ; filevarattdef(fout,"time_bounds" ,time_bounds) ; fout->time_bounds=(/time_bounds/) ; if( isfilepresent("domain_xy.nc") )then fxy = addfile("domain_xy.nc","r") fout->x = fxy->x fout->y = fxy->y fout->ProjectionCoordinateSystem = fxy->ProjectionCoordinateSystem proj4=fxy->ProjectionCoordinateSystem@proj4 esri_pe_string=fxy->ProjectionCoordinateSystem@esri_pe_string else print("domain_xy.nc file is missing") end if ; v1d = (/ndtooned((/f1->T2D(0,:,:)/))/) if(indexFlag .eq. "false") then findexOut = addfile("index_out.nc","c") index = ind(ismissing(v1d) ) findexOut->index = index system("mv index_out.nc index.nc") exit end if v1dgfs = (/ndtooned((/f2->T2D(0,:,:)/))/) v1d(index) = v1dgfs(index) delete(v1dgfs) T2D = onedtond(v1d,(/1,dimlat,dimlon/)) delete(v1d) if(any(ismissing(T2D)))then T2D = where(ismissing(T2D),f2->T2D,T2D) print("WARNING: combine.ncl, Missing value is found in T2D, replaced with RAP data.") ; ; Email SPA ; system("echo -e '\nWARNING: "+hrrrFile+" missing value is found in T2D, replaced with RAP data.' >> $mailMsgNCL") ; end if T2D@_FillValue = floattointeger(t2d_missing) T2D@units="K" T2D@long_name="2-m Air Temperature" T2D@standard_name = "air_temperature" ; T2D@scale_factor = t2d_scale_fd; T2D@add_offset = t2d_offsetd; T2D@valid_range = t2d_valid_range ; T2D@missing_value = floattointeger(t2d_missing) ; T2D@proj4=proj4 T2D@grid_mapping=grid_mapping T2D@esri_pe_string=esri_pe_string T2D@remap=remap ; T2D@cell_methods=cell_methods T2D!0="time" T2D!1="y" T2D!2="x" filevardef(fout, "T2D" ,typeof(T2D) ,getvardims(T2D)) filevarattdef(fout,"T2D" ,T2D) fout->T2D = (/T2D/) delete(T2D) ; Downward long wave radiation using GFS. v1d =(/ ndtooned(f2->LWDOWN) /) if((.not.isatt(f2->LWDOWN,"scale_factor")).and.(.not.isatt(f2->LWDOWN,"add_offset")))then v1d1=floattointeger(getDataScaled(v1d,lwdown_scale_f,lwdown_offset)) else v1d1=v1d end if delete(v1d) LWDOWN = onedtond(v1d1,(/1,dimlat,dimlon/)) delete(v1d1) if(any(ismissing(LWDOWN)))then LWDOWN = where(ismissing(LWDOWN),100.,LWDOWN) print("WARNING: combine.ncl, Missing value is found in LWDOWN, using RAP data.") ; ; Email SPA ; system("echo -e '\nWARNING: "+gfsFile+" missing value is found in LWDOWN, using RAP data.' >> $mailMsgNCL") ; end if LWDOWN@_FillValue=floattointeger(lwdown_missing) LWDOWN@units="W m-2" LWDOWN@long_name="Surface downward long-wave radiation flux" LWDOWN@standard_name = "surface_downward_longwave_flux" ; LWDOWN@scale_factor = lwdown_scale_fd; LWDOWN@add_offset = lwdown_offsetd ; LWDOWN@valid_range = lwdown_valid_range ; LWDOWN@missing_value = floattointeger(lwdown_missing) ; LWDOWN@proj4=proj4 LWDOWN@grid_mapping=grid_mapping LWDOWN@esri_pe_string=esri_pe_string LWDOWN@remap=remap ; LWDOWN@cell_methods=cell_methods LWDOWN!0="time" LWDOWN!1="y" LWDOWN!2="x" filevardef(fout, "LWDOWN" ,typeof(LWDOWN) ,getvardims(LWDOWN)) filevarattdef(fout,"LWDOWN" ,LWDOWN) fout->LWDOWN = (/LWDOWN/) delete(LWDOWN) v1d = (/ndtooned((/f1->Q2D(0,:,:)/)) /) v1dgfs = (/ndtooned((/f2->Q2D(0,:,:)/))/) v1d(index) = v1dgfs(index) delete(v1dgfs) Q2D = onedtond(v1d,(/1,dimlat,dimlon/)) delete(v1d) if(any(ismissing(Q2D)))then Q2D = where(ismissing(Q2D),f2->Q2D,Q2D) print("WARNING: combine.ncl, Missing value is found in Q2D, using RAP data.") ; ; Email SPA ; system("echo -e '\nWARNING: "+hrrrFile+" missing value is found in Q2D, using RAP data.' >> $mailMsgNCL") ; end if Q2D@_FillValue=floattointeger(q2d_missing) Q2D@units="kg kg-1" Q2D@long_name="2-m Specific Humidity, dimensionless ratio of the mass of water vapor (kg) to the total mass of the system (kg)" Q2D@standard_name = "surface_specific_humidity" ; Q2D@scale_factor = q2d_scale_fd; Q2D@add_offset = q2d_offsetd ; Q2D@valid_range = q2d_valid_range ; Q2D@missing_value = floattointeger(q2d_missing) ; Q2D@proj4=proj4 Q2D@grid_mapping=grid_mapping Q2D@esri_pe_string=esri_pe_string Q2D@remap=remap ; Q2D@cell_methods=cell_methods Q2D!0="time" Q2D!1="y" Q2D!2="x" Q2D!0="time" Q2D!1="y" Q2D!2="x" filevardef(fout, "Q2D" ,typeof(Q2D) ,getvardims(Q2D)) filevarattdef(fout,"Q2D" ,Q2D) fout->Q2D = (/Q2D/) delete(Q2D) v1d = (/ ndtooned(f1->U2D) /) v2dgfs = (/f2->U2D/) v1dgfs = ndtooned(v2dgfs) v1d(index) = v1dgfs(index) delete(v1dgfs) v1d1=floattointeger(getDataScaled(v1d,u2d_scale_f,u2d_offset)) delete(v1d) U2D = onedtond(v1d1,(/1,dimlat,dimlon/)) delete(v1d1) if(any(ismissing(U2D)))then U2D = where(ismissing(U2D),v2dgfs,U2D) print("WARNING: combine.ncl, Missing value is found in U2D, using RAP data.") ; ; Email SPA ; system("echo -e '\nWARNING: "+hrrrFile+" missing value is found in U2D, using RAP data.' >> $mailMsgNCL") ; end if delete(v2dgfs) U2D@_FillValue=floattointeger(u2d_missing) U2D@units="m s-1" U2D@long_name="10-m U-component of wind" U2D@standard_name = "x_wind" ; U2D@scale_factor = u2d_scale_fd; U2D@add_offset = u2d_offsetd ; U2D@valid_range = u2d_valid_range ; U2D@missing_value = floattointeger(u2d_missing) ; U2D@proj4=proj4 U2D@grid_mapping=grid_mapping U2D@esri_pe_string=esri_pe_string U2D@remap=remap ; U2D@cell_methods=cell_methods U2D!0="time" U2D!1="y" U2D!2="x" filevardef(fout, "U2D" ,typeof(U2D) ,getvardims(U2D)) filevarattdef(fout,"U2D" ,U2D) fout->U2D = (/U2D/) delete(U2D) v1d = (/ ndtooned(f1->V2D) /) v2dgfs = (/f2->V2D/) v1dgfs = ndtooned(v2dgfs) v1d(index) = v1dgfs(index) delete(v1dgfs) v1d1=floattointeger(getDataScaled(v1d,v2d_scale_f,v2d_offset)) delete(v1d) V2D = onedtond(v1d1,(/1,dimlat,dimlon/)) delete(v1d1) if(any(ismissing(V2D)))then V2D = where(ismissing(V2D),v2dgfs,V2D) print("WARNING: combine.ncl, Missing value is found in V2D, using RAP data.") ; ; Email SPA ; system("echo -e '\nWARNING: "+hrrrFile+" missing value is found in V2D, using RAP data.' >> $mailMsgNCL") ; end if delete(v2dgfs) V2D@_FillValue=floattointeger(v2d_missing) V2D@units="m s-1" V2D@long_name="10m V-component of wind" V2D@standard_name = "y_wind" ; V2D@scale_factor = v2d_scale_fd; V2D@add_offset = v2d_offsetd ; V2D@valid_range = v2d_valid_range ; V2D@missing_value = floattointeger(v2d_missing) ; V2D@proj4=proj4 V2D@grid_mapping=grid_mapping V2D@esri_pe_string=esri_pe_string V2D@remap=remap ; V2D@cell_methods=cell_methods V2D!0="time" V2D!1="y" V2D!2="x" filevardef(fout, "V2D" ,typeof(V2D) ,getvardims(V2D)) filevarattdef(fout,"V2D" ,V2D) fout->V2D = (/V2D/) delete(V2D) v1d = (/ ndtooned((/f1->PSFC(0,:,:)/)) /) v1dgfs = (/ ndtooned((/f2->PSFC(0,:,:)/)) /) v1d(index) = v1dgfs(index) delete(v1dgfs) PSFC = onedtond(v1d,(/1,dimlat,dimlon/)) delete(v1d) if(any(ismissing(PSFC)))then PSFC = where(ismissing(PSFC),f2->PSFC,PSFC) print("WARNING: combine.ncl, Missing value is found in PSFC, using RAP data.") ; ; Email SPA ; system("echo -e '\nWARNING: "+hrrrFile+" missing value is found in PSFC, using RAP data.' >> $mailMsgNCL") ; end if PSFC@_FillValue=floattointeger(psfc_missing) PSFC@units="Pa" PSFC@long_name = "Surface Pressure" PSFC@standard_name = "air_pressure" ; PSFC@scale_factor = psfc_scale_fd; PSFC@add_offset = psfc_offsetd; PSFC@valid_range = psfc_valid_range ; PSFC@missing_value = floattointeger(psfc_missing) ; PSFC@proj4=proj4 PSFC@grid_mapping=grid_mapping PSFC@esri_pe_string=esri_pe_string PSFC@remap=remap ; PSFC@cell_methods=cell_methods PSFC!0="time" PSFC!1="y" PSFC!2="x" filevardef(fout, "PSFC" ,typeof(PSFC) ,getvardims(PSFC)) filevarattdef(fout,"PSFC" ,PSFC) fout->PSFC = (/PSFC/) delete(PSFC) v1d = (/ ndtooned(f1->RAINRATE) /) v2dgfs = (/f2->RAINRATE/) v1dgfs = ndtooned(v2dgfs) v1d(index) = v1dgfs(index) delete(v1dgfs) RAINRATE = onedtond(v1d,(/1,dimlat,dimlon/)) delete(v1d) if(any(ismissing(RAINRATE)))then RAINRATE = where(ismissing(RAINRATE),v2dgfs,RAINRATE) print("WARNING: combine.ncl, Missing value is found in RAINRATE, using RAP data.") ; ; Email SPA ; system("echo -e '\nWARNING: "+hrrrFile+" missing value is found in RAINRATE, using RAP data.' >> $mailMsgNCL") ; end if delete(v2dgfs) RAINRATE@_FillValue=floattointeger(rainrate_missing) RAINRATE@units="mm s-1" RAINRATE@long_name="Surface Precipitation Rate" RAINRATE@standard_name = "precipitation_flux" ; RAINRATE@valid_range = rainrate_valid_range ; RAINRATE@missing_value = rainrate_missing ; RAINRATE@proj4=proj4 RAINRATE@grid_mapping=grid_mapping RAINRATE@esri_pe_string=esri_pe_string RAINRATE@remap=remap ; RAINRATE@cell_methods=cell_methods_rain RAINRATE!0="time" RAINRATE!1="y" RAINRATE!2="x" filevardef(fout, "RAINRATE" ,typeof(RAINRATE) ,getvardims(RAINRATE)) filevarattdef(fout,"RAINRATE" ,RAINRATE) fout->RAINRATE = (/RAINRATE/) delete(RAINRATE) v1d = (/ndtooned((/f1->SWDOWN(0,:,:)/))/) v1dgfs = (/ndtooned((/f2->SWDOWN(0,:,:)/))/) v1d(index) = v1dgfs(index) delete(v1dgfs) SWDOWN = onedtond(v1d,(/1,dimlat,dimlon/)) delete(v1d) if(any(ismissing(SWDOWN)))then SWDOWN = where(ismissing(SWDOWN),f2->SWDOWN,SWDOWN) print("WARNING: combine.ncl, Missing value is found in SWDOWN, using RAP data.") ; ; Email SPA ; system("echo -e '\nWARNING: "+hrrrFile+" missing value is found in SWDOWN, using RAP data.' >> $mailMsgNCL") ; end if SWDOWN@_FillValue=floattointeger(swdown_missing) SWDOWN@units="W m-2" SWDOWN@long_name="Surface downward short-wave radiation flux" SWDOWN@standard_name = "surface_downward_shortwave_flux" ; SWDOWN@scale_factor = swdown_scale_fd; SWDOWN@add_offset = swdown_offsetd; SWDOWN@valid_range = swdown_valid_range ; SWDOWN@missing_value = floattointeger(swdown_missing) ; SWDOWN@proj4=proj4 SWDOWN@grid_mapping=grid_mapping SWDOWN@esri_pe_string=esri_pe_string SWDOWN@remap=remap ; SWDOWN@cell_methods=cell_methods SWDOWN!0="time" SWDOWN!1="y" SWDOWN!2="x" filevardef(fout, "SWDOWN" ,typeof(SWDOWN) ,getvardims(SWDOWN)) filevarattdef(fout,"SWDOWN" ,SWDOWN) fout->SWDOWN = (/SWDOWN/) delete(SWDOWN) end