;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ; 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 "../ncl-beta/contributed_beta.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl" ; 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/shea_util.ncl" begin ; smooth 9 point smooth = False ;smooth = True ;smooth_twice = False smooth_twice = True ; Define some variables that we will use in the formatting of input and output file names ; Optional, specific to my directory structure Dnum="01" Dsize="12" case="20141125" year="2014" month="11" day="25" starthour="18" model="RAP" PBL="MYNN2" micro="THOM" dom="NE" ; filename="wrfout_d03_2016-01-23_00:00:00.nc" ; filename="wrfout_d03_2016-01-23_00:00:00.nc" ; filename="wrfout_d02_2016-01-22_18:00:00_mar_con.nc" filename="/data1/wrf/02_18_2013/wrfout_d04_2013-02-08_00:00:00_mp8" ; filename="/data1/wrf/02_18_2013/wrfout_d03_2013-02-08_00:00:00_mp8" ; filename="/data1/wrf/02_18_2013/wrfout_d02_2013-02-08_00:00:00_mp8" ; filename="wrfout_d01_2016-01-22_18:00:00_1dom.nc" ; The WRF ARW input file. ; This needs to have a ".nc" appended, so just do it. ; a = addfile("../Research/"+case+"_Case/"+model+"_"+PBL+"_"+micro+"/wrfout_d"+Dnum+"_"+year+"-"+month+"-"+day+"_"+starthour+":00:00.nc","r") a = addfile(filename,"r") ; We generate plots, but what kind do we prefer? ; type = "x11" ;type = "png" type = "pdf" ; type = "ps" ; type = "ncgm" ; Set some basic resources res = True ; res@MainTitle = "12 km WRF: "+model+"-"+PBL+"-"+micro res@MainTitle = "1.3 km WRF" ; res@MainTitle = "4.0 km WRF" ; res@MainTitle = "12.0 km WRF" lats = (/ 38.00, 42.50 /) lons = (/ -76.0, -68.75 /) loc = wrf_user_ll_to_ij(a, lons, lats, True) ; loc(0,;) is west-east (x) ; loc(1,:) is south-north (y) ; subtract one since we want to use it as an index in NCL x_start = loc(0,0) - 1 x_end = loc(0,1) - 1 y_start = loc(1,0) - 1 y_end = loc(1,1) - 1 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 mpres@ZoomIn = True mpres@Xstart = x_start mpres@Ystart = y_start mpres@Xend = x_end mpres@Yend = y_end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; What times and how many time steps are in the data set? ; find ntimes times = wrf_user_list_times(a) ; get times in the file ntimes = dimsizes(times) ; number of times in the file print ("ntimes = " + ntimes) ; wks = gsn_open_wks(type,"front_sbm_1.3km") wks = gsn_open_wks(type,"fronto_mp8_02_08_2013_1.3km") ; wks = gsn_open_wks(type,"fronto_mp8_02_08_2013_4.0km") ; wks = gsn_open_wks(type,"fronto_mp8_02_08_2013_12.0km") ; wks = gsn_open_wks(type,"front_sbm_4.0km") ; do it = 6,6,1 ; do it = 3,7,1 ; do it = 12,36,1 do it = 46,54,1 ; do it = 47,47,1 gsn_define_colormap(wks,"MPL_YlOrRd") ; 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 lat = wrf_user_getvar(a,"lat",it) lon = wrf_user_getvar(a,"lon",it) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; 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) z_dec@units="dam" z_dec@short_name="Height" ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Calculate desired function that involves gradients ; barry xlat = wrf_user_getvar(a,"XLAT",it) xlong = wrf_user_getvar(a,"XLONG",it) src_lat = a->XLAT(0,:,:) ;;---Change (maybe) src_lon = a->XLONG(0,:,:) ;;---Change (maybe) printVarSummary(xlat) ; ; 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. ;If you are regridding to a rectilinear grid, then set: ;opt@DstGridType = "rectilinear" ; If you are regridding to an unstructured grid, set: ; opt@DstGridType = "unstructured" Opt@DstGridType = "rectilinear" 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) ;SrcFileName is just the name you want to give the intermediate NetCDF file that gets created in order ;to describe your source grid to the ESMF_RegridWeightGen program. If you don't set this, the ;intermediate file will be called "source_grid_file.nc". ;I suggest not setting this at all, or, if you do set it, call it something like "WRF_SCRIP_file.nc". ; Opt@SrcInputFileName = src_file ; optional, but good idea ; Opt@SrcInputFileName = a ; 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 if (smooth_twice)then do i_smooth=1,9 th_plane=smth9(th_plane, 0.50, 0.25, False) ; Heavy local smoothing u_plane= smth9(u_plane, 0.50, .25, False) ; light local smoothing v_plane=smth9(v_plane, 0.50, .25, False) ; light local smoothing z_plane= smth9(z_plane, 0.50, .25, False) ; light local smoothing end do end if if (smooth)then th_plane=smth9(th_plane, 0.50, 0.25, False) ; Heavy local smoothing u_plane= smth9(u_plane, 0.50, .25, False) ; light local smoothing v_plane=smth9(v_plane, 0.50, .25, False) ; light local smoothing z_plane= smth9(z_plane, 0.50, .25, False) ; light local smoothing end if 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 ; barry del_th = grad_latlon_cfd(th_regrid,dst_lat,dst_lon,False,False) dTdy = del_th[0] ; latitude gradient of var_regrid dTdx = del_th[1] ; longitude gradient dTdx@units = "K/m" dTdy@units = "K/m" ; barry del_u = grad_latlon_cfd(u_regrid,dst_lat,dst_lon,False,False) dudy = del_u[0] ; latitude gradient of var_regrid dudx = del_u[1] ; longitude gradient del_v = grad_latlon_cfd(v_regrid,dst_lat,dst_lon,False,False) dvdy = del_th[0] ; latitude gradient of var_regrid dvdx = del_th[1] ; longitude gradient ; Do the calculation with those gradient variables mag_del_theta = sqrt(dTdx^2 + dTdy^2) ; Magnitude of del theta ;rh90 = where(rh_plane.gt.90,rh_plane,rh_plane@_FillValue) ; only want to shade RH values greater than 90% ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Other formulation of FGEN ; Copied from GEMPAK appendix B1, page 6 conversion = 10^5 * 10800 ; Scale to SPC Mesoanalysis units of K/100 km/3 hr divergence = dudx + dvdy shearing = dudy + dvdy stretching = dudx - dvdy delta = 0.5 * atan(shearing/stretching) beta = asin((-cos(delta) * dTdx - sin(delta) * dTdy) / mag_del_theta) resultant_def = sqrt(shearing^2 + stretching^2) pett_fgen = 0.5 * conversion * mag_del_theta * (resultant_def * cos(2. * beta) - divergence) wrf_smooth_2d( pett_fgen, 3 ) pett_fgen@units = "K/100 km/3 hr" pett_fgen@short_name = "Petterssen Frontogenesis" ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Set some plotting resources for our final map ; Plotting options for fgen opts = res ; opts@units="K/100 km/3 hr" ; opts@UnitLabel = "Frontogenesis" ; opts@units="K/100 km/3 hr" opts@FieldTitle = "Frontogenesis" opts@UnitLabel = "K/100 km/3 hr" opts@cnLevelSelectionMode = "ExplicitLevels" opts@cnFillOn = True ; opts@cnLevels = (/ 0, 1,2., \ ; 3, 4.,5.,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20/) ; 3, 4.,5.,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20/) opts@cnFillColors = (/"White","White","AntiqueWhite","AntiqueWhite1","AntiqueWhite2","AntiqueWhite3",\ "olivedrab","olivedrab1","olivedrab2", \ "olivedrab3","SpringGreen","SpringGreen2","SpringGreen4", \ "Yellow","Yellow2","Yellow4","Orange","Red","Red3","HotPink3","HotPink1","HotPink","Violet"/) opts@cnLevels = (/ 0, 1,3., \ 5, 7.,9.,11,13,15,17,19,21,23,25,30,35,40,45,50,75,100/) opts@gsnAddCyclic = False opts@cnFillMode = "RasterFill" ; for large data sets var_zoom = pett_fgen(y_start:y_end,x_start:x_end) ; contour_fgen = wrf_contour(a,wks,pett_fgen,opts) contour_fgen = wrf_contour(a,wks,var_zoom,opts) delete(opts) ; Plotting options for rh optsrh = res optsrh@cnLineColor = "green3" optsrh@cnMonoFillColor = True optsrh@cnFillColor = "green3" optsrh@cnFillOn = True optsrh@cnLinesOn = False optsrh@cnLineLabelsOn = False optsrh@cnConstFLabelPerimOn = False ;optsrh@ContourParameters = (/ 90., 100., 10./) optsrh@cnFillOpacityF = 0.2 ; 20% opacity optsrh@cnLineLabelBackgroundColor = -1 ;optsrh@gsnContourLineThicknessesScale = 1.5 optsrh@cnLabelMasking = True optsrh@cnInfoLabelOn = False optsrh@cnLineLabelPerimOn = False optsrh@lbLabelBarOn = False ; Do not display label bar optsrh@gsnAddCyclic = False ; var_zoom = rh(y_start:y_end,x_start:x_end) ;contour_rh = wrf_contour(a,wks,rh90,optsrh) ;contour_rh = wrf_contour(a,wks,var_zoom,optsrh) delete(optsrh) ; Plotting options for Z opts = res opts@FieldTitle = "Heights" opts@UnitLabel = "dam" opts@cnLineColor = "Black" opts@cnHighLabelsOn = False ;opts@cnHighLabelFontColor = "blue" opts@cnLowLabelsOn = False ;opts@cnLowLabelFontColor = "red" opts@cnConstFLabelPerimOn = False opts@ContourParameters = (/ 210., 390., 3. /) opts@cnLineLabelBackgroundColor = -1 opts@gsnContourLineThicknessesScale = 3 opts@cnLabelMasking = True opts@cnInfoLabelOn = False opts@cnLineLabelPerimOn = False opts@gsnAddCyclic = False var_zoom = z_plane(y_start:y_end,x_start:x_end) ; contour_z = wrf_contour(a,wks,z_plane,opts) contour_z = wrf_contour(a,wks,var_zoom,opts) delete(opts) ; Plotting options for SLP opts = res opts@cnLineColor = "Black" opts@cnHighLabelsOn = True opts@cnHighLabelFontColor = "blue" opts@cnLowLabelsOn = True opts@cnLowLabelFontColor = "red" opts@cnConstFLabelPerimOn = False opts@ContourParameters = (/ 900., 1100., 4. /) opts@cnLineLabelBackgroundColor = -1 ; opts@gsnContourLineThicknessesScale = 3.0 opts@cnLabelMasking = True opts@cnInfoLabelOn = False opts@cnLineLabelPerimOn = False opts@gsnAddCyclic = False var_zoom = slp(y_start:y_end,x_start:x_end) ; contour_psl = wrf_contour(a,wks,slp,opts) ; contour_psl = wrf_contour(a,wks,var_zoom,opts) delete(opts) ;---Write variables to file for Mary ; system("rm -f file_for_mary.nc") ; fout = addfile("file_for_mary.nc","c") ; fout->pett_fgen = pett_fgen ; fout->z_plane = z_plane ; fout->xlat = xlat ; fout->xlong = xlong ; MAKE PLOTS plot = wrf_map_overlays(a,wks,(/contour_fgen,contour_z/),pltres,mpres) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; end do ; END OF TIME LOOP end