; load "$USHnwm/scale/set_scale_factor.ncl" ; filename1="$USHnwm/scale/set_scale_factor.ncl" ; if(isfilepresent(filename1)) then ; loadscript(filename1) ; scaleFlag="True" ; else ; print("Error: scale factor file doesn't exist: " + filename1) ; scaleFlag="False" ; end if begin w1=WEIGHT1/3.0 file1=getenv("FILE_in1") file2=getenv("FILE_in2") file3=getenv("FILE_in3") VALIDTIME=getenv("VALIDTIME") w2=(3-WEIGHT1)/3.0 ff1=addfile(file1,"r") ff2=addfile(file2,"r") if(isfilepresent(file3)) then system("/bin/rm -f "+file3) end if ncdf=addfile(file3,"c") globalAtt = True globalAtt@model_output_valid_time = VALIDTIME globalAtt@model_initialization_time = ff1@model_initialization_time globalAtt@NWM_version_number=nwm_version globalAtt@model_output_type="forcing" globalAtt@model_configuration="medium_range" globalAtt@model_total_valid_times=240 fileattdef( ncdf, globalAtt ) T2D0 = ff1->T2D dd = dimsizes(T2D0) ndd=dimsizes(dd) print("dd="+dd+" ndd="+ndd) dimlat = dd(ndd-2) dimlon = dd(ndd-1) delete(T2D0) ; dimNames = (/"time", "y", "x","reference_time","nv"/) ; dimSizes = (/ 1 , dimlat, dimlon , 1, 2/) ; dimUnlim = (/ False , False, False,False,False/) dimNames = (/"time", "y", "x","reference_time"/) dimSizes = (/ 1 , dimlat, dimlon , 1/) dimUnlim = (/ False , False, False,False/) filedimdef(ncdf,dimNames,dimSizes,dimUnlim) time=tointeger(w1*(ff1->time)+w2*(ff2->time)) time!0="time" time@long_name="valid output time" ; time@bounds="time_bounds" time@standard_name="time" time@calendar=ff1->time@calendar time@units=ff1->time@units time@valid_max=ff1->time@valid_max time@valid_min=ff1->time@valid_min filevardef(ncdf, "time" ,typeof(time) ,getvardims(time)) filevarattdef(ncdf,"time" ,time) ncdf->time=(/time/) ; if( isfilevar(ff1,"reference_time") ) then ncdf->reference_time=ff1->reference_time end if ; time_bounds=new((/1,2/),integer) time_bounds!0="time" time_bounds!1="nv" time_bounds@units=ff1->time_bounds@units time_bounds(0,0)=(/ff1->reference_time(0)/) time_bounds(0,1)=doubletoint(time(0)) delete(time_bounds@_FillValue) ; filevardef(ncdf, "time_bounds" ,typeof(time_bounds) ,getvardims(time_bounds)) ; filevarattdef(ncdf,"time_bounds" ,time_bounds) ; ncdf->time_bounds=(/time_bounds/) ; if( isfilepresent("domain_xy.nc") )then fxy = addfile("domain_xy.nc","r") ncdf->x = fxy->x ncdf->y = fxy->y ncdf->ProjectionCoordinateSystem = fxy->ProjectionCoordinateSystem esri_pe_string=fxy->ProjectionCoordinateSystem@esri_pe_string else print("domain_xy.nc file is missing") end if T2D0 = new((/1,dimlat,dimlon/),float) temp1=ff1->T2D temp2=ff2->T2D if ((isatt(temp1,"scale_factor")).and.(isatt(temp2,"scale_factor"))) T2D = new((/1,dimlat,dimlon/),integer) T2D0 = w1*(getDataUnscaled(temp1,doubletofloat(temp1@scale_factor),doubletofloat(temp1@add_offset))) \ +w2*(getDataUnscaled(temp2,doubletofloat(temp2@scale_factor),doubletofloat(temp2@add_offset))) T2D=floattointeger(getDataScaled(T2D0,t2d_scale_f,t2d_offset)) T2D@scale_factor = t2d_scale_fd; T2D@add_offset = t2d_offsetd ; T2D@valid_range = t2d_valid_range ; T2D@missing_value = floattointeger(t2d_missing) ; T2D@_FillValue = floattointeger(t2d_missing) else T2D = new((/1,dimlat,dimlon/),typeof(ff1->T2D)) T2D = w1*(temp1)+w2*(temp2) T2D@missing_value=getVarFillValue(temp1) T2D@_FillValue=getVarFillValue(temp1) end if T2D@units="K" T2D@long_name="2-m Air Temperature" T2D@standard_name = "air_temperature" ; 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(ncdf, "T2D" ,typeof(T2D) ,getvardims(T2D)) filevarattdef(ncdf,"T2D" ,T2D) ncdf->T2D = (/T2D/) delete(T2D) delete(T2D0) delete(temp1) delete(temp2) Q2D0 = new((/1,dimlat,dimlon/),float) temp1=ff1->Q2D temp2=ff2->Q2D if ((isatt(temp1,"scale_factor")).and.(isatt(temp2,"scale_factor"))) Q2D = new((/1,dimlat,dimlon/),integer) Q2D0 = w1*(getDataUnscaled(temp1,doubletofloat(temp1@scale_factor),doubletofloat(temp1@add_offset))) \ +w2*(getDataUnscaled(temp2,doubletofloat(temp2@scale_factor),doubletofloat(temp2@add_offset))) Q2D=floattointeger(getDataScaled(Q2D0,q2d_scale_f,q2d_offset)) Q2D@scale_factor = q2d_scale_fd; Q2D@add_offset = q2d_offsetd ; Q2D@valid_range = q2d_valid_range ; Q2D@missing_value = floattointeger(q2d_missing) ; Q2D@_FillValue = floattointeger(q2d_missing) else Q2D = new((/1,dimlat,dimlon/),typeof(temp1)) Q2D = w1*temp1+w2*temp2 Q2D@missing_value=temp1@missing_value Q2D@_FillValue=temp1@_FillValue end if Q2D@units="kg kg-1" Q2D@long_name="2-m Specific Humidity, the dimensionless ratio of the mass of water vapor (kg) to the total mass of the system (kg)" Q2D@standard_name = "surface_specific_humidity" ; 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" filevardef(ncdf, "Q2D" ,typeof(Q2D) ,getvardims(Q2D)) filevarattdef(ncdf,"Q2D" ,Q2D) ncdf->Q2D = (/Q2D/) delete(Q2D) delete(Q2D0) delete(temp1) delete(temp2) U2D0 = new((/1,dimlat,dimlon/),float) temp1=ff1->U2D temp2=ff2->U2D if ((isatt(temp1,"scale_factor")).and.(isatt(temp2,"scale_factor"))) U2D = new((/1,dimlat,dimlon/),integer) U2D0 = w1*(getDataUnscaled(temp1,doubletofloat(temp1@scale_factor),doubletofloat(temp1@add_offset))) \ +w2*(getDataUnscaled(temp2,doubletofloat(temp2@scale_factor),doubletofloat(temp2@add_offset))) U2D=floattointeger(getDataScaled(U2D0,u2d_scale_f,u2d_offset)) U2D@scale_factor = u2d_scale_fd; U2D@add_offset = u2d_offsetd ; U2D@valid_range = u2d_valid_range ; U2D@missing_value = floattointeger(u2d_missing) ; U2D@_FillValue = floattointeger(u2d_missing) else U2D = new((/1,dimlat,dimlon/),typeof(ff1->U2D)) U2D = w1*(ff1->U2D)+w2*(ff2->U2D) U2D@missing_value=temp1@missing_value U2D@_FillValue=temp1@_FillValue end if U2D@units="m s-1" U2D@long_name="10-m U-component of wind" U2D@standard_name = "x_wind" ; 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(ncdf, "U2D" ,typeof(U2D) ,getvardims(U2D)) filevarattdef(ncdf,"U2D" ,U2D) ncdf->U2D = (/U2D/) delete(U2D) delete(U2D0) delete(temp1) delete(temp2) V2D0 = new((/1,dimlat,dimlon/),float) temp1=ff1->V2D temp2=ff2->V2D if ((isatt(temp1,"scale_factor")).and.(isatt(temp2,"scale_factor"))) V2D = new((/1,dimlat,dimlon/),integer) V2D0 = w1*(getDataUnscaled(temp1,doubletofloat(temp1@scale_factor),doubletofloat(temp1@add_offset))) \ +w2*(getDataUnscaled(temp2,doubletofloat(temp2@scale_factor),doubletofloat(temp2@add_offset))) V2D=floattointeger(getDataScaled(V2D0,v2d_scale_f,v2d_offset)) V2D@scale_factor = v2d_scale_fd; V2D@add_offset = v2d_offsetd ; V2D@valid_range = v2d_valid_range ; V2D@missing_value = floattointeger(v2d_missing) ; V2D@_FillValue = floattointeger(v2d_missing) else V2D = new((/1,dimlat,dimlon/),typeof(ff1->V2D)) V2D = w1*(ff1->V2D)+w2*(ff2->V2D) V2D@missing_value=temp1@missing_value V2D@_FillValue=temp1@_FillValue end if V2D@units="m s-1" V2D@long_name="10-m V-component of wind" V2D@standard_name = "y_wind" ; 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(ncdf, "V2D" ,typeof(V2D) ,getvardims(V2D)) filevarattdef(ncdf,"V2D" ,V2D) ncdf->V2D = (/V2D/) delete(V2D) delete(V2D0) delete(temp1) delete(temp2) PSFC0 = new((/1,dimlat,dimlon/),float) temp1=ff1->PSFC temp2=ff2->PSFC if ((isatt(temp1,"scale_factor")).and.(isatt(temp2,"scale_factor"))) PSFC = new((/1,dimlat,dimlon/),integer) PSFC0 = w1*(getDataUnscaled(temp1,doubletofloat(temp1@scale_factor),doubletofloat(temp1@add_offset))) \ +w2*(getDataUnscaled(temp2,doubletofloat(temp2@scale_factor),doubletofloat(temp2@add_offset))) PSFC=floattointeger(getDataScaled(PSFC0,psfc_scale_f,psfc_offset)) PSFC@scale_factor = psfc_scale_fd; PSFC@add_offset = psfc_offsetd ; PSFC@valid_range = psfc_valid_range ; PSFC@missing_value = floattointeger(psfc_missing) ; PSFC@_FillValue = floattointeger(psfc_missing) else PSFC = new((/1,dimlat,dimlon/),typeof(temp1)) PSFC = w1*temp1+w2*temp2 PSFC@missing_value=temp1@missing_value PSFC@_FillValue=temp1@_FillValue end if PSFC@units="Pa" PSFC@long_name = "Surface Pressure" PSFC@standard_name = "air_pressure" ; 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(ncdf, "PSFC" ,typeof(PSFC) ,getvardims(PSFC)) filevarattdef(ncdf,"PSFC" ,PSFC) ncdf->PSFC = (/PSFC/) delete(PSFC) delete(PSFC0) delete(temp1) delete(temp2) LWDOWN0 = new((/1,dimlat,dimlon/),float) temp1=ff1->LWDOWN temp2=ff2->LWDOWN if ((isatt(temp1,"scale_factor")).and.(isatt(temp2,"scale_factor"))) LWDOWN = new((/1,dimlat,dimlon/),integer) LWDOWN0 = w1*(getDataUnscaled(temp1,doubletofloat(temp1@scale_factor),doubletofloat(temp1@add_offset))) \ +w2*(getDataUnscaled(temp2,doubletofloat(temp2@scale_factor),doubletofloat(temp2@add_offset))) LWDOWN=floattointeger(getDataScaled(LWDOWN0,lwdown_scale_f,lwdown_offset)) LWDOWN@scale_factor = lwdown_scale_fd; LWDOWN@add_offset = lwdown_offsetd ; LWDOWN@valid_range = lwdown_valid_range ; LWDOWN@missing_value = floattointeger(lwdown_missing) ; LWDOWN@_FillValue = floattointeger(lwdown_missing) else LWDOWN = new((/1,dimlat,dimlon/),typeof(ff1->LWDOWN)) LWDOWN = w1*temp1+w2*temp2 LWDOWN@missing_value=temp1@missing_value LWDOWN@_FillValue=temp1@_FillValue end if LWDOWN@units="W m-2" LWDOWN@long_name="Surface downward long-wave radiation flux" LWDOWN@standard_name = "surface_downward_longwave_flux" ; 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(ncdf, "LWDOWN" ,typeof(LWDOWN) ,getvardims(LWDOWN)) filevarattdef(ncdf,"LWDOWN" ,LWDOWN) ncdf->LWDOWN = (/LWDOWN/) delete(LWDOWN) delete(LWDOWN0) delete(temp1) delete(temp2) SWDOWN0 = new((/1,dimlat,dimlon/),float) temp1=ff1->SWDOWN temp2=ff2->SWDOWN if ((isatt(temp1,"scale_factor")).and.(isatt(temp2,"scale_factor"))) SWDOWN = new((/1,dimlat,dimlon/),integer) SWDOWN0 = w1*(getDataUnscaled(temp1,doubletofloat(temp1@scale_factor),doubletofloat(temp1@add_offset))) \ +w2*(getDataUnscaled(temp2,doubletofloat(temp2@scale_factor),doubletofloat(temp2@add_offset))) SWDOWN=floattointeger(getDataScaled(SWDOWN0,swdown_scale_f,swdown_offset)) SWDOWN@scale_factor = swdown_scale_fd; SWDOWN@add_offset = swdown_offsetd ; SWDOWN@valid_range = swdown_valid_range ; SWDOWN@missing_value = floattointeger(swdown_missing) ; SWDOWN@_FillValue = floattointeger(swdown_missing) else SWDOWN = new((/1,dimlat,dimlon/),typeof(temp1)) SWDOWN = w1*temp1+w2*temp2 SWDOWN@missing_value=temp1@missing_value SWDOWN@_FillValue=temp1@_FillValue end if SWDOWN@units="W m-2" SWDOWN@long_name="Surface downward short-wave radiation flux" SWDOWN@standard_name = "surface_downward_shortwave_flux" ; 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(ncdf, "SWDOWN" ,typeof(SWDOWN) ,getvardims(SWDOWN)) filevarattdef(ncdf,"SWDOWN" ,SWDOWN) ncdf->SWDOWN = (/SWDOWN/) delete(SWDOWN) delete(SWDOWN0) delete(temp1) delete(temp2) RAINRATE = new((/1,dimlat,dimlon/),float) RAINRATE = w1*(ff1->RAINRATE)+w2*(ff2->RAINRATE) RAINRATE@valid_range = rainrate_valid_range ; RAINRATE@missing_value = rainrate_missing ; RAINRATE@_FillValue = rainrate_missing ; RAINRATE@description = "RAINRATE" RAINRATE@long_name = "Surface Precipitation Rate" RAINRATE@standard_name = "precipitation_flux" ; RAINRATE@units = "mm s-1" 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(ncdf, "RAINRATE" ,typeof(RAINRATE) ,getvardims(RAINRATE)) filevarattdef(ncdf,"RAINRATE" ,RAINRATE) ncdf->RAINRATE = (/RAINRATE/) delete(RAINRATE) end