; --------------------------------------------------------------------------- ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Program Name: NAM2Hydro_forcing_ESMFregrid.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 perform regridding. ;
; ;
; Input: nam file, weighting function, output file ;
; ;
; Output: regridded file ;
; ;
; For non-fatal errors output is witten to $DATA/logs ;
; ;
; Author(s)/Contact(s): Linlin Pan, lpan@ucar.edu ;
; Origination Sept., 2017 ;
; ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;----------------------------------------------------------------------
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/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl"
;load "./ESMF_regridding.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"
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
;----------------------------------------------------------------------
; Source, destination , and weight filenames for generating
; ESMF_regridding weights
;----------------------------------------------------------------------
wgtFileName_conserve = wgtFileName_in
if ( .not.isfilepresent( dstGridName ) ) then
print( " ... source grid file not found : "+ dstGridName )
exit
end if
;----------------------------------------------------------------------
; read in source and destination grid netCDF files
;----------------------------------------------------------------------
;--- destination grid data
dstfile = addfile( dstGridName ,"r")
dlon3d=dstfile->XLONG_M ;(USER'S NEED TO CONFIRM THIS VARIABLE IS WHAT IS IN THEIR DATA FILE)
dlat3d=dstfile->XLAT_M ;(USER'S NEED TO CONFIRM THIS VARIABLE IS WHAT IS IN THEIR DATA FILE)
dlon2d=dlon3d(0,:,:)
dlat2d=dlat3d(0,:,:)
dims=dimsizes(dlat2d)
outFile = getenv("outFile")
;
;dg NEED TO EDIT INPUT DATA TIME INTERVAL HERE...
;
dt=3600.0 ;forcing data timestep in seconds... (USER'S MUST ENTER/CONFIRM THIS IS SET TO DATA TIMESTEP (or INTERVAL))
flag=0 ;WRF - flag for removing accum precip... (DO NOT CHANGE THIS)
;----------------------------------------------------------------------
; Open source data files to be regridded...
;----------------------------------------------------------------------
;
outdir = "./output_files/" ; directory where output forcing data will be placed. set to dirm for overwriting the original file
if(.not. isfilepresent(outdir)) then
system("mkdir "+outdir)
end if
if(isfilepresent(outdir+"/"+outFile)) then
system("rm "+outdir+"/"+outFile)
end if
srcfilename = getenv ("srcFile")
datfils = systemfunc ("/bin/ls -1 "+srcfilename) ;list of file names
num_datfils = dimsizes(datfils)
init_time = getenv("init_time") ; outFile data
valid_time = getenv("valid_time") ; outFile data
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)
print("yyww check valid_time :"+valid_time)
print("yyww check init_time :"+init_time)
wgtFileName = wgtFileName_in
opt = True
opt@WgtFileName = wgtFileName
opt@CopyVarAtts = True
opt@CopyVarCoords = False
opt@Debug = True
do ifil = 0,num_datfils-1,1 ; loop through datafiles one at a time
datfile = addfile( datfils(ifil), "r")
print( " ... Open input file : "+ datfils(ifil) )
names = getfilevarnames(datfile)
;----------------------------------------------------------------------
; Temporary output
;----------------------------------------------------------------------
ncdf= addfile(outdir+"/"+outFile,"c")
system("rm -f test.nc")
; ncdf->lat = dlat2d ;output lat
; ncdf->lon = dlon2d ;output lon
globalAtt = True
globalAtt@NWM_version_number=nwm_version
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"
globalAtt@model_output_type="forcing"
globalAtt@model_configuration="short_range"
globalAtt@model_total_valid_times=60
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))
yeari=tointeger(init_dtime(0))
monthi=tointeger(init_dtime(1))
dayi=tointeger(init_dtime(2))
houri=tointeger(init_dtime(3))
houre=houri+60
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))
time = tointeger(cd_inv_calendar(year,month,day,hour,0,0,units, 0))
fileattdef( ncdf, globalAtt )
;----------------------------------------------------------------------
; Processing...no further mods should be required...
;----------------------------------------------------------------------
;do v=6,6
; print("lv_HTGL2 = "+ datfile->lv_HTGL2)
; print("lv_HTGL5 = "+ datfile->lv_HTGL5)
; print("lv_HTGL9 = "+ datfile->lv_HTGL9)
T2D0 = ESMF_regrid_with_weights(datfile->TMP_P0_L103_GME0(0,:,:) , wgtFileName, opt)
dd = dimsizes(T2D0)
print("dd="+dd)
dimlat = dd(0)
dimlon = dd(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(ncdf,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@calendar = "standard" ;
time@valid_max=time_max
time@valid_min=time_min
time@units = units ;
filevardef(ncdf, "time" ,typeof(time) ,getvardims(time))
filevarattdef(ncdf,"time" ,time)
ncdf->time=(/time/)
;
yeari=tointeger(init_dtime(0))
monthi=tointeger(init_dtime(1))
dayi=tointeger(init_dtime(2))
houri=tointeger(init_dtime(3))
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
;
filevardef(ncdf, "reference_time" ,typeof(reference_time) ,getvardims(reference_time))
filevarattdef(ncdf,"reference_time" ,reference_time)
ncdf->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)=toint(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
proj4=fxy->ProjectionCoordinateSystem@proj4
else
print("domain_xy.nc file is missing")
end if
T2D = new((/1,dimlat,dimlon/),float)
T2D(0,:,:)=T2D0
delete(T2D0)
T2D@missing_value=-1.e+33
T2D@_FillValue=9.96921e+36
T2D@units="K"
T2D@long_name="2-m Air Temperature"
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)
Q2D = new((/1,dimlat,dimlon/),float)
Q2D(0,:,:) = ESMF_regrid_with_weights(datfile->SPFH_P0_L103_GME0(0,:,:) , wgtFileName, opt)
Q2D@missing_value=-1.e+33
Q2D@_FillValue=9.96921e+36
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!0="time"
Q2D!1="y"
Q2D!2="x"
filevardef(ncdf, "Q2D" ,typeof(Q2D) ,getvardims(Q2D))
filevarattdef(ncdf,"Q2D" ,Q2D)
ncdf->Q2D = (/Q2D/)
delete(Q2D)
;10 meter wind
U2D0 = new((/1,dimlat,dimlon/),float)
U2D0(0,:,:) = ESMF_regrid_with_weights(datfile->UGRD_P0_L103_GME0(0,:,:) , wgtFileName, opt)
U2D=floattointeger(getDataScaled(U2D0,u2d_scale_f,u2d_offset))
U2D@_FillValue=floattointeger(u2d_missing)
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@units="m s-1"
U2D@long_name="10-m U-component of 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)
V2D0 = new((/1,dimlat,dimlon/),float)
V2D0(0,:,:) = ESMF_regrid_with_weights(datfile->VGRD_P0_L103_GME0(0,:,:) , wgtFileName, opt)
V2D=floattointeger(getDataScaled(V2D0,v2d_scale_f,v2d_offset))
V2D@_FillValue=floattointeger(v2d_missing)
V2D@long_name="10-m 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@units="m s-1"
V2D@long_name="10-m V-component of 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)
PSFC = new((/1,dimlat,dimlon/),float)
PSFC(0,:,:) = ESMF_regrid_with_weights(datfile->PRES_P0_L1_GME0 , wgtFileName, opt)
PSFC@missing_value=-1.e+33
PSFC@_FillValue=9.96921e+36
PSFC@units="Pa"
PSFC@long_name = "Surface Pressure"
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)
RAINRATE = new((/1,dimlat,dimlon/),float)
RAINRATE = 0.0
RAINRATE@_FillValue=(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@description = "RAINRATE"
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)
do i = 0, dimsizes( names ) - 1
if(names(i) .eq. "PRATE_P0_L1_GME0") then
ytmp = ESMF_regrid_with_weights(datfile->PRATE_P0_L1_GME0, wgtFileName, opt)
RAINRATE(0,:,:) =(/ytmp/)
delete(ytmp)
end if
if(names(i) .eq. "PRATE_P8_L1_GME0_avg") then
ytmp = ESMF_regrid_with_weights(datfile->PRATE_P8_L1_GME0_avg, wgtFileName, opt)
RAINRATE(0,:,:) =(/ytmp/)
delete(ytmp)
end if
if(names(i) .eq. "PRATE_P8_L1_GME0_avg6h") then
ytmp = ESMF_regrid_with_weights(datfile->PRATE_P8_L1_GME0_avg6h, wgtFileName, opt)
RAINRATE(0,:,:) =(/ytmp/)
delete(ytmp)
end if
if(names(i) .eq. "PRATE_P8_L1_GME0_avg3h") then
ytmp = ESMF_regrid_with_weights(datfile->PRATE_P8_L1_GME0_avg3h, wgtFileName, opt)
RAINRATE(0,:,:) =(/ytmp/)
delete(ytmp)
end if
end do
RAINRATE = where(RAINRATE<0,0,RAINRATE)
ncdf->RAINRATE = (/RAINRATE/)
delete(RAINRATE)
LWDOWN0 = new((/1,dimlat,dimlon/),float)
do i = 0, dimsizes( names ) - 1
if(names(i) .eq. "DLWRF_P0_L1_GME0") then
LWDOWN0(0,:,:) = ESMF_regrid_with_weights(datfile->DLWRF_P0_L1_GME0, wgtFileName, opt)
LWDOWN=floattointeger(getDataScaled(LWDOWN0,lwdown_scale_f,lwdown_offset))
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(ncdf, "LWDOWN" ,typeof(LWDOWN) ,getvardims(LWDOWN))
filevarattdef(ncdf,"LWDOWN" ,LWDOWN)
ncdf->LWDOWN = (/LWDOWN/)
end if
if(names(i) .eq. "DLWRF_P8_L1_GME0_avg") then
LWDOWN0(0,:,:) = ESMF_regrid_with_weights(datfile->DLWRF_P8_L1_GME0_avg, wgtFileName, opt)
LWDOWN=floattointeger(getDataScaled(LWDOWN0,lwdown_scale_f,lwdown_offset))
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(ncdf, "LWDOWN" ,typeof(LWDOWN) ,getvardims(LWDOWN))
filevarattdef(ncdf,"LWDOWN" ,LWDOWN)
ncdf->LWDOWN = (/LWDOWN/)
end if
end do
delete(LWDOWN)
delete(LWDOWN0)
SWDOWN = new((/1,dimlat,dimlon/),float)
do i = 0, dimsizes( names ) - 1
if(names(i) .eq. "DSWRF_P0_L1_GME0") then
SWDOWN(0,:,:) = ESMF_regrid_with_weights(datfile->DSWRF_P0_L1_GME0 , wgtFileName, opt)
SWDOWN@missing_value=-1.e+33
SWDOWN@_FillValue=9.96921e+36
SWDOWN@units="W m-2"
SWDOWN@long_name="Surface downward short-wave radiation flux"
SWDOWN!0="time"
SWDOWN!1="y"
SWDOWN!2="x"
filevardef(ncdf, "SWDOWN" ,typeof(SWDOWN) ,getvardims(SWDOWN))
filevarattdef(ncdf,"SWDOWN" ,SWDOWN)
ncdf->SWDOWN = (/SWDOWN/)
end if
end do
delete(SWDOWN)
; WEASD = ESMF_regrid_with_weights(datfile->WEASD_P0_L1_GME0, wgtFileName, opt)
; ncdf->WEASD = WEASD
; delete(WEASD)
end do ; end do for file loop
end