;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ; This script computes temperature advection from a wrfout .nc file. ; ; It can be modified to compute anything that requires a gradient. ; ; Computing the horizontal gradient of a variable on a regional grid, ; such as WRF, requires that the user: ; ; 1. First convert the WRF grid from its native curvilinear ; format to rectilinear format using ESMF_regrid for every ; variable which will need to have its gradient calculated, ; or act on a gradient variable. ; This is built in to this script for T in Kelvin, and u and v ; wind components. ; (For more information on ESMF_regrid, see: ; http://www.ncl.ucar.edu/Document/Functions/ESMF/ESMF_regrid.shtml) ; ; 2. Either have NCL version 6.4 installed, or download the newest ; contributed.ncl (http://www.cgd.ucar.edu/~shea/contributed.ncl_beta_640) ; and load it instead of the older contributed.ncl ; This is necessary because the computation of gradients will call ; on the function grad_latlon_cfd, which is only available in this ; latest version. ; ; 3. Call the grad_latlon_cfd function with the new rectilinear variable ; grids as arguments and a variable as the output assignment. It will ; return the y- and x-gradient as the 0th and 1st element of the array, ; respectively. ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" load "local_bin/contributed.ncl_beta_640" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl" begin ; Define some variables that we will use in the formatting of input and output file names ; Optional, specific to my directory structure filename="wrfout_d01_2016-01-22_18:00:00.nc" ; The WRF ARW input file. ; This needs to have a ".nc" appended, so just do it. a = addfile(filename,"r") ; We generate plots, but what kind do we prefer? ; type = "x11" type = "png" ; type = "ps" ; type = "ncgm" ; Set some basic resources res = True res@MainTitle = "12 km WRF 300 hPa PV, 700 hPa Pett. Fgen" pltres = True mpres = True mpres@mpGeophysicalLineColor = "Black" mpres@mpNationalLineColor = "Black" mpres@mpUSStateLineColor = "Black" mpres@mpNationalLineColor = "Black" mpres@mpGridLineColor = 0 mpres@mpLimbLineColor = "Black" mpres@mpPerimLineColor = "Black" mpres@mpUSStateLineThicknessF = 1.0 mpres@mpGeophysicalLineThicknessF = 1.0 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; What times and how many time steps are in the data set? times = wrf_user_getvar(a,"times",-1) ; get all times in the file ntimes = dimsizes(times) ; number of times in the file do it = 12,ntimes-1,12 ; TIME LOOP wks = gsn_open_wks(type,"wrf_processed_maps") gsn_define_colormap(wks,"MPL_Purples") ; Overwrite the standard color map print("Working on time: " + times(it) ) res@TimeLabel = times(it) ; Set Valid time to use on plots ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; First get the variables we will need slp = wrf_user_getvar(a,"slp",it) ; slp wrf_smooth_2d( slp, 3 ) ; smooth slp th = wrf_user_getvar(a,"th",it) ; 3D theta ; td = wrf_user_getvar(a,"td",it) ; 3D td u = wrf_user_getvar(a,"ua",it) ; 3D U at mass points v = wrf_user_getvar(a,"va",it) ; 3D V at mass points p = wrf_user_getvar(a, "pressure",it) ; pressure is our vertical coordinate rh = wrf_user_getvar(a, "rh",it) ; relative humidity z = wrf_user_getvar(a,"z",it) ; Full model height in meters z_dec = z/10. ; Height in decameters pv = wrf_user_getvar(a,"pvo",it) ; Potential vorticity ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Interpolate to isobaric planes pressure = 700. th_plane = wrf_user_intrp3d(th,p,"h",pressure,0.,False) z_plane = wrf_user_intrp3d(z_dec,p,"h",pressure,0.,False) rh_plane = wrf_user_intrp3d(rh,p,"h",pressure,0.,False) u_plane = wrf_user_intrp3d( u,p,"h",pressure,0.,False) v_plane = wrf_user_intrp3d( v,p,"h",pressure,0.,False) pv_plane = wrf_user_intrp3d(pv,p,"h",300,0.,False) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Regrid variable(s) to compute gradients on the regional grid ; Only the variables that will have gradients computed need to be regridded ;---WRF file containing source grid src_file="wrfout_d01_2016-01-22_18:00:00_new.nc" ; src_file = a sfile = addfile(src_file,"r") ;---Get variable to regrid ;var = tk_plane ;;---Change (likely) src_lat = sfile->XLAT(0,:,:) ;;---Change (maybe) src_lon = sfile->XLONG(0,:,:) ;;---Change (maybe) ; ; Create the destination lat/lon rectilinear grid. ; ; Note: the code below is just an example of how you might ; create a rectilinear grid to regrid to. The code below ; is using the source grid to determine the min/max ; lat/lon values, and the number of points, and then ; generating 1D lat/lon arrays based on this. ; ; It is possible that dst_lat and dst_lon could be read ; off of another file. ; src_dims = dimsizes(src_lat) nlat = src_dims(0) nlon = src_dims(1) dst_lat = fspan(min(src_lat),max(src_lat),nlat) ;;---Change (likely) dst_lon = fspan(min(src_lon),max(src_lon),nlon) ;;---Change (likely) ;---Set up regridding options Opt = True ;---"bilinear" is the default. "patch" and "conserve" are other options. Opt@InterpMethod = "bilinear" ;;---Change (maybe) Opt@WgtFileName = "wrfout_temp_rect.nc" Opt@SrcGridLat = src_lat ; source grid Opt@SrcGridLon = src_lon Opt@SrcRegional = True ;;--Change (maybe) Opt@SrcInputFileName = src_file ; optional, but good idea Opt@DstGridLat = dst_lat ; destination grid Opt@DstGridLon = dst_lon Opt@DstRegional = True ;;--Change (maybe) Opt@ForceOverwrite = True Opt@PrintTimings = False Opt@Debug = True th_regrid = ESMF_regrid(th_plane,Opt) ; Do the regridding u_regrid = ESMF_regrid(u_plane,Opt) ; Do the regridding v_regrid = ESMF_regrid(v_plane,Opt) ; Do the regridding ; printVarSummary(tk_regrid) ; printVarSummary(u_regrid) ; printVarSummary(v_regrid) end do ; END OF TIME LOOP end