;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; wrfout_to_cf.ncl ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; -NCL script to read an ARW wrfout NetCDF file on staggered model ; grid and to output unstaggered values in NetCDF CF compliant format. ; -Extensive help in maintaining and adding to the script was provide by ; Matt Higgins. ; ; command syntax: ; ncl 'file_in="wrfout.nc"' 'file_out="wrfpost.nc"' wrfout_to_cf.ncl ; ; -The NCL script is executed by the above command syntax. Alternatively, ; the file_out and file_in can be set in the script and there is then no ; need to specify it at the command prompt. ; -The values that are to be included in the output are determined by ; setting several attribute variables to True (include) or False (skip). ; These attribute variables are set at around line 160 in this script. ; -Setting the overall variable of a group of variables to false will exclude ; all of the variables in that group from the output file. For example, ; setting the variable out2dMet to False means that all 2dMet variables will ; not be included. If out2dMet is set to true, all of the individual variable ; attributes (i.e. out2dMet@T_sfc, out2dMet@T_2m, etc.) that are also set to ; true will be included in the output file. ; ; Support information: ; This script is semi-supported on a time available basis. If ; you come across an error, or have an idea for an improvement, please send ; an email. I will update the script on a time when I have time. Send all ; inquiries or questions to: Mark Seefeldt - mark.seefeldt@colorado.edu ; ; Release Notes: v2.0.0 ; -release notes for all versions can be found at: ; http://foehn.colorado.edu/wrfout_to_cf/ ; -v2.0.0 ; NOTE: A significant change was made to re-name all mixing ratio ; variables to 'r_' from 'q'. Examples: ; Before: q_2m Now: r_v_2m ; q_e r_v_e ; q_p r_v_p ; q_cloud r_cloud ; q_rain_p r_rain_p ; This change was made with the addition of specific humidty as a ; moisture output option. The decision was made to go with the more ; community accepted variable of 'r' to represent mixing ratio. ; The previous values of q_2m, q_e, and q_p now represent specific ; humidity. ; ; -Added specific humidity as an output at 2m, eta levels, and pressure levels. ; -All measures of moisture (Td, RH, r_v, q) are available at 2m, eta levels, ; and pressure levels. ; -Added winds rotated to true earth coord. at 2m, eta, and pressure levels. ; -Added wind direction (Earth) and wind speed at 2m, eta, and pressure levels. ; -Added potential vorticity and absolute vorticity at eta and pressure levels. ; -Changed several variable names to more simples forms. ; ex: SHFlx and LHFlx to SH and LH. ; -Added a CLWRF section of a selection of CLWRF variables. ; -Reorganized the code for a simpler and easy to follow layout. ; ; TODO: In an upcoming release there will be an option to use CMIP variable ; names in the output file. There are placeholders for this option ; that are included below, but this is not an active feature at this ; point. ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; load in the libraries load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; procedure to process the attributes in creating CF compliant WRF output procedure assignVarAttCoord(x:numeric, time[*]:numeric, vert[*]:numeric, \ fl_vert:numeric) ; x:numeric -variable to process the attributes ; time[*]:numeric -array with time values for adding coordinates ; vert[*]:numeric -array with vertical values for adding coordinates ; Note: set to 0 if no vertical coordinate ; fl_vert:numeric -flag indicating vertical coordinate type ; 0 = no vertical coordinate (x,y only) ; 1 = pressure ; 2 = eta ; 3 = soil ; MissingValue -assigned missing value attribute begin ; assign the default missing value MissingValue = 1e20 ; set time for all variables x!0 = "time" x&time = time ; set the vertical coordinate depending on fl_vert if (fl_vert .eq. 1) then ;pressure as vertical coordinate x!1 = "pressure" x&pressure = vert ;x@missing_value = MissingValue x@_FillValue = MissingValue end if if (fl_vert .eq. 2) then ;eta as vertical coordinate x!1 = "eta" x&eta = vert x@_FillValue = MissingValue end if if (fl_vert .eq. 3) then ;soil as vertical coordinate x!1 = "soil" x&soil = vert x@_FillValue = MissingValue end if ; set the horizontal coordinates if (fl_vert .eq. 0) then ;no vertical coordinate x!1 = "south_north" x!2 = "west_east" ;x@missing_value = MissingValue x@_FillValue = MissingValue else ;with vertical coordinate x!2 = "south_north" x!3 = "west_east" ;x@missing_value = MissingValue x@_FillValue = MissingValue end if ; set the mapping coordinates x@coordinates = "lon lat" end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; start the primary wrfout_to_cf.ncl program begin ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; configuration settings ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; set the units for time TimeUnits = "hours since 2001-01-01 00:00:00" ; set the values for pressure to be interpolated pressure = (/1000.,850.,700.,500.,300./) ; set the limits for the output range ; 0 = beginning of dataset ; 9999 = end of dataset ; Note: remember that the array is zero-based limTime = (/0,9999/) limS_N = (/0,9999/) limW_E = (/0,9999/) limPres = (/0,9999/) limEta = (/0,9999/) limSoil = (/0,9999/) ; set default values for file_in, dir_in, and file_out, if not specified if (.not.isvar("file_in")) then file_in = "wrfout_d02_2005-11-19_00:00:00" end if if (.not.isvar("dir_in")) then dir_in = "./" end if if (.not.isvar("file_out")) then file_out = "wrfpost.nc" end if if (.not.isvar("dir_out")) then dir_out = "./" end if ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;; set the flags for selecting variables to be included ;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Note: Not all of the indicated values below can be extracted / ; converted from a given wrfout file. Some of the fields are ; included with only specific WRF physics options. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; output settings axis = True ;one-dimensional coordinate fields projection = False ;CF projection info with fields outPtop = False ;include Ptop in the output file ;outCMIP = False ;use CMIP variable names in output ;MissingValue = -999999 ;missing value ; Note: MissingValue is currently assigned in the procedure assignVarAttCoord ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; set the netcdf file global attributes fileAtt = True fileAtt@creation_date = systemfunc("date") fileAtt@institution = "University of Colorado at Boulder - CIRES" fileAtt@created_by = "Mark Seefeldt - mark.seefeldt@colorado.edu" fileAtt@notes = "Created with NCL script: wrfout_to_cf.ncl v2.0" fileAtt@source = file_in fileAtt@Conventions = "CF 1.6, Standard Name Table v19" fileAtt@title = file_out ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; time / date variables outDateTime = True ;include a yyyymmddhh field outUTCDate = True ;include yr,mo,dy,hr,mn fields ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; two-dimensional near-surface / surface met variables out2dMet = True out2dMet@T_sfc = False ;temperature at the surface out2dMet@slp = True ;sea-level pressure - using WRF-NCL out2dMet@u_10m_tr = False ;u wind - rotated to earth - at 10m out2dMet@v_10m_tr = False ;v wind - rotated to earth - at 10m out2dMet@precip_g = False ;total grid scale precipitation ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; three-dimensional upper-level (eta levels) meteorology variables outEta = False ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; three-dimensional upper-level (pressure levels) meteorology variables outPressure = False ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; two-dimensional surface energy budget / radiation variables out2dRadFlx = False out2dRadFlx@SH = False ;SH flux - upward - sfc - instant out2dRadFlx@LH = False ;LH flux - upward - sfc - instant ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; three-dimensional soil variables outSoil = False ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; two-dimensional CLimate WRF (CLWRF) variables ; (http://www.meteo.unican.es/en/software/clwrf) ; -The CLWRF variables are activated by compiler flags within configiure.wrf. ; -The variables are outputted to auxiliary history output file(s). outCLWRF = False ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; open the input netcdf file (wrfout file) wrfout = addfile(dir_in+file_in+".nc","r") ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; time coordinate ; -the time in wrfout is in an odd character format TimeChar = wrfout->Times ; -determine the number of dimensions for time DimTimeChar = dimsizes(TimeChar) nTime = DimTimeChar(0) ; -convert the wrfout time to a CF compliant time ; "hours since 1901-01-01 00:00:00" time_in = wrf_times_c(TimeChar, 1) ; -create an array indicating the year, month, day, hour, minute, second utc_date = floattoint(ut_calendar(time_in, 0)) ; -create the final variable for time with the units selected time = (/ut_inv_calendar(utc_date(:,0), utc_date(:,1), utc_date(:,2), \ utc_date(:,3), utc_date(:,4), utc_date(:,5), \ TimeUnits, 0)/) ;time time@long_name = "Time" time@standard_name = "time" time@units = TimeUnits time@calendar = "standard" time!0 = "time" time&time = time utc_date!0 = "time" ;utc_date utc_date&time = time year = utc_date(:,0) year@long_name = "Year" year!0 = "time" year&time = time month = utc_date(:,1) month@long_name = "Month" month!0 = "time" month&time = time day = utc_date(:,2) day@long_name = "Day" day!0 = "time" day&time = time hour = utc_date(:,3) hour@long_name = "Hour" hour!0 = "time" hour&time = time minute = utc_date(:,4) minute@long_name = "Minutes" minute!0 = "time" minute&time = time ; -convert the wrfout time to a DateTime integer for easy reading if (outDateTime) then DateTime = (/wrf_times_c(TimeChar, 3)/) ;time DateTime@long_name = "Date and Time" DateTime!0 = "time" DateTime&time = time end if ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; vertical variables / coordinates ; Note: pressure levels are assigned in the beginning section if (outPressure) then nPressure = dimsizes(pressure) ;pressure vertical coordinate pressure@long_name = "Pressure Levels" pressure@standard_name = "air_pressure" pressure@units = "hPa" pressure@positive = "down" pressure!0 = "pressure" pressure&pressure = pressure end if if (outEta .or. outPressure) eta = (/wrfout->ZNU(0,:)/) ;eta values on half-levels (mass) nEta = dimsizes(eta) eta@long_name = "Eta Levels (mass points)" eta@standard_name = "atmosphere_sigma_coordinate" eta@units = "1" eta@positive = "down" eta@formula_terms = "sigma: eta ps: p_sfc ptop: p_top" eta!0 = "eta" eta&eta = eta end if if (outSoil) then soil = (/wrfout->ZS(0,:)/) ;depths of center of soil layers nSoil = dimsizes(soil) soil@long_name = "Soil Levels (depth)" soil@standard_name = "depth" soil@units = "m" soil@positive = "down" soil!0 = "soil" soil&soil = soil end if ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; one-dimensional general model variables if (outPtop .or. outEta) then p_top_in = (/wrfout->P_TOP/)/100. ;pressure at top of model ;in some files P_TOP has two dimensions, in some it has one dimension if ((dimsizes(dimsizes(p_top_in))) .eq. 2) then p_top = p_top_in(0,0) else p_top = p_top_in(0) end if p_top@long_name = "Pressure at Top of the Model" p_top@standard_name = "air_pressure" p_top@units = "hPa" end if ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; two-dimensional mapping variables lat = (/wrfout->XLAT(0,:,:)/) ;lat (mass) DimLat = dimsizes(lat) nS_N = DimLat(0) ;S_N dimension nW_E = DimLat(1) ;W_E dimension lat@long_name = "Latitude" lat@standard_name = "latitude" lat@units = "degrees_north" lat!0 = "south_north" lat!1 = "west_east" lon = (/wrfout->XLONG(0,:,:)/) ;lon (mass) lon@long_name = "Longitude" lon@standard_name = "longitude" lon@units = "degrees_east" lon!0 = "south_north" lon!1 = "west_east" Z_sfc = (/wrfout->HGT(0,:,:)/) ;Z_sfc Z_sfc@long_name = "Terrain Height" Z_sfc@standard_name = "height" Z_sfc@units = "m" Z_sfc@coordinates = "lon lat" Z_sfc!0 = "south_north" Z_sfc!1 = "west_east" ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; one-dimensional coordinate system if (axis) then south_north = ispan(0,nS_N-1,1) ;south_north south_north@long_name = "y-coordinate in Cartesian system" south_north@axis = "Y" south_north@units = "m" south_north!0 = "south_north" west_east = ispan(0,nW_E-1,1) ;west_east west_east@long_name = "x-coordinate in Cartesian system" west_east@axis = "X" west_east@units = "m" west_east!0 = "west_east" end if ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;check the limits for the output arrays, set 9999 to end of dataset if (limTime(1) .eq. 9999) then limTime(1) = nTime-1 end if if (limS_N(1) .eq. 9999) then limS_N(1) = nS_N-1 end if if (limW_E(1) .eq. 9999) then limW_E(1) = nW_E-1 end if if (outPressure) then if (limPres(1) .eq. 9999) then limPres(1) = nPressure-1 end if end if if (outEta) then if (limEta(1) .eq. 9999) then limEta(1) = nEta-1 end if end if if (outSoil) then if (limSoil(1) .eq. 9999) then limSoil(1) = nSoil-1 end if end if ;create filename and open post-processed netCDF file if isfilepresent(dir_out+file_out) then system ("rm "+dir_out+file_out) ;remove any pre-exisiting file end if wrfpost = addfile(dir_out+file_out,"c") ;create new netCDF file filedimdef (wrfpost, "time", nTime, True) ; establish a variable for a new line in the attributes nl = integertochar(10) ; newline character ; create the global attributes fileattdef(wrfpost, fileAtt) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;; write post-processed WRF data to netCDF file ;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ; check to see if the output file is to follow CMIP guidelines ; if (outCMIP .ne. True) then ;use standard wrfout_to_cf variable names in the wrfpost output file ; -date and time variables wrfpost->time=time(limTime(0):limTime(1)) if (outDateTime) then wrfpost->DateTime=DateTime(limTime(0):limTime(1)) end if if (outUTCDate) then wrfpost->year = year(limTime(0):limTime(1)) wrfpost->month = month(limTime(0):limTime(1)) wrfpost->day = day(limTime(0):limTime(1)) wrfpost->hour = hour(limTime(0):limTime(1)) wrfpost->minute = minute(limTime(0):limTime(1)) end if ; -vertical coordinate variables if (outPressure) then wrfpost->pressure=pressure(limPres(0):limPres(1)) end if if (outEta) then wrfpost->eta=eta(limEta(0):limEta(1)) end if if (outSoil) then wrfpost->soil=soil(limSoil(0):limSoil(1)) end if ; -one-dimensional coordinate variables ; This inclusion of this section of output results in errors when using ; NCL 6.0 and later. Not certain as to why. ; if (axis) then ; wrfpost->south_north=south_north(limS_N(0):limS_N(1)) ; wrfpost->west_east=west_east(limW_E(0):limW_E(1)) ; end if ; -one-dimensional general model variables if (outPtop .or. outEta) then wrfpost->p_top=p_top end if ; -two-dimensional mapping variables wrfpost->lat=lat(limS_N(0):limS_N(1),limW_E(0):limW_E(1)) wrfpost->lon=lon(limS_N(0):limS_N(1),limW_E(0):limW_E(1)) wrfpost->Z_sfc=Z_sfc(limS_N(0):limS_N(1),limW_E(0):limW_E(1)) ; -two-dimensional near-surface / surface met variables ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; two-dimensional near-surface / surface met variables ; -retrieved directly from the wrfout file - or ; -derived/diagnostic fields using wrf_user_getvar if (out2dMet) then if (out2dMet@T_sfc) then T_sfc = (/wrfout->TSK/) ;T_sfc T_sfc@long_name = "Temperature at the Surface" T_sfc@standard_name = "surface_temperature" T_sfc@units = "K" assignVarAttCoord(T_sfc,time,0,0) wrfpost->T_sfc=T_sfc delete(T_sfc) end if if (out2dMet@slp) then ;slp slp = new( (/ nTime, nS_N, nW_E /), float) slp@long_name = "Sea-Level Pressure" slp@standard_name = "air_pressure_at_sea_level" slp@units = "hPa" assignVarAttCoord(slp,time,0,0) print(dimsizes(slp)) do n = 0, 3 print(n) ;slp = (/wrf_user_getvar(wrfout,"slp",n)/) ;print(dimsizes((/wrf_user_getvar(wrfout,"slp",n)/))) slp(n,:,:)=(/wrf_user_getvar(wrfout,"slp",n)/) ;wrfpost->slp(n,:,:)=(/wrf_user_getvar(wrfout,"slp",n)/) end do wrfpost->slp =slp ;wrfpost->slp=slp delete(slp) end if ; if (out2dMet@u_10m_tr .or. out2dMet@v_10m_tr) then ; U_10m = new( (/ 2, nTime, nS_N, nW_E /), float) ; do n = 0, nTime-1 ; print(n) ; U_10m = (/wrf_user_getvar(wrfout,"uvmet10",n)/) ; end do ; u_10m_tr = (/U_10m(0,:,:,:)/) ;u_10m_tr ; u_10m_tr@long_name = "u-Component of wind at 10 m (Earth)" ; u_10m_tr@standard_name = "eastward_wind" ; u_10m_tr@units = "m s-1" ; assignVarAttCoord(u_10m_tr,time,0,0) ; wrfpost->u_10m_tr=u_10m_tr ; v_10m_tr = (/U_10m(1,:,:,:)/) ;v_10m_tr ; v_10m_tr@long_name = "v-Component of wind at 10 m (Earth)" ; v_10m_tr@standard_name = "northward_wind" ; v_10m_tr@units = "m s-1" ; assignVarAttCoord(v_10m_tr,time,0,0) ; wrfpost->v_10m_tr=v_10m_tr ; delete([/U_10m,u_10m_tr,v_10m_tr/]) ; end if if (out2dMet@precip_g) then precip_g = (/wrfout->RAINNC/) ;precip_g precip_g@long_name = "Accumulated Total Grid Scale Precipitation" precip_g@standard_name = "large_scale_precipitation_amount" precip_g@units = "mm" assignVarAttCoord(precip_g,time,0,0) wrfpost->precip_g=precip_g delete(precip_g) end if end if ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; two-dimensional surface energy budget / radiation variables if (out2dRadFlx) then if (out2dRadFlx@SH) then SH = (/wrfout->HFX/) ;SH flux - upward - sfc - instant SH@long_name = "Sensible Heat Flux - Upward - at surface - instant" SH@standard_name = "surface_upward_sensible_heat_flux" SH@units = "W m-2" assignVarAttCoord(SH,time,0,0) wrfpost->SH=SH delete(SH) end if if (out2dRadFlx@LH) then LH = (/wrfout->LH/) ;LH flux - upward - sfc - instant LH@long_name = "Latent Heat Flux - Upward - at surface - instant" LH@standard_name = "surface_upward_latent_heat_flux" LH@units = "W m-2" assignVarAttCoord(LH,time,0,0) wrfpost->LH=LH delete(LH) end if end if ; end if ; ; check to see if the output file is to follow CMIP guidelines ; if (outCMIP .eq. True) then ; ; end if end