;====================================================================== ; ESMF_regrid_37.ncl ; ; Concepts illustrated: ; - Extract needed variables from a WRF file [curvilinear grid] ; - Create a rectilinear grid with 'nominally' the same curvilinear resolution ; - Use ESMF [bilinear] to regrid the WRF curvilinear u and v ; to the 'nominal' rectilinear grid ; - Use the weight file generated when regridding 'u' for the 'v' variable ; - Use NCL's ' uv2dv_cfd ' to compute the divergence on the rectilinear grid ; - Use ESMF to regrid the ivergence on the nominal rectilinear grid ; back onto the original WRF curvilinear grid ; - Plot ;====================================================================== ; Note: Technically, u and v must be regridded as a pair because there ; is some reorientation when wind components are interpolated ; to different latitude and longitude locations.. ;====================================================================== ;---No need to load these after 6.1.2 ;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/wrf/WRFUserARW.ncl" ;---Must load the 'ESMF-regridding.ncl' library thru 6.4.0 ;---Automatically loaded from 6.5.0 onward load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl" ;====================================================================== begin ;--- Addfiles a = addfile("/media/lara/MyPassport/SIMULACIONES/STC16_parametrizaciones/STC16/nuevoD02/output_20141020/wrfout_d02_2014-10-19_12___00___00_411", "r") ;;;;Change Date;;;; x = 1530 ;;time = 3 time = 0 ; one time step on this file ;--- Get variables lat2d = wrf_user_getvar(a, "XLAT", time) ; (south_north,west_east) lon2d = wrf_user_getvar(a, "XLONG",time) u = wrf_user_getvar(a,"ua",time) ; (bottom_top,south_north,west_east) v = wrf_user_getvar(a,"va",time) printVarSummary(u) printMinMax(u,0) print("----") printVarSummary(v) printMinMax(v,0) print("----") vo = wrf_user_getvar(a,"avo",time) ; calculate avo for all times Ft = a->F F= Ft(2:36,:,:) wr = vo-F ; relative vorticity printVarSummary(wr) printMinMax(wr,0) print("----") ;---Create the destination rectilinear lat[*]/lon[*] arrays ;---Nominally, the same resolution dims = dimsizes(lat2d) nlat = dims(0) nlon = dims(1) lat = fspan(min(lat2d), max(lat2d) ,nlat) ; approximately same resolution lon = fspan(min(lon2d), max(lon2d) ,nlon) lat!0 = "lat" lon!0 = "lon" ;exit() ;================================================================================ ; Part A: WRF (curvilinear) to Rectilinear regrid variable and generate weights ;================================================================================ ;Interpolate the WRF u, v and relative vorticity to a rectilinear grid with as fine [or finer] resoultion as the WRF grid ;---Set up options for regridding InterpMethod= "bilinear" ; define interpolation method ;---Title and filename options Opt = True Opt@SrcTitle = "WRF grid" ; optional Opt@WgtFileName = "WRF_to_Rect.WgtFile_"+InterpMethod+"."+nlat+"_x_"+nlon+".nc" ;---Source grid options Opt@SrcFileName = "WRF.SCRIP_grid_description."+nlat+"_x_"+nlon+".nc" ; Name of source and Opt@SrcRegional = True Opt@SrcGridLat = lat2d Opt@SrcGridLon = lon2d ;---Destination grid options Opt@DstFileName = "Rectilinear.SCRIP_grid_description."+nlat+"_x_"+nlon+".nc" ; destination files Opt@DstGridType = "rectilinear" Opt@DstGridLat = lat Opt@DstGridLon = lon Opt@DstRegional = True ;---Specify other options Opt@ForceOverwrite = True Opt@InterpMethod = InterpMethod ;---Perform the regrid: WRF ==> rectilinear (_reclin) u_reclin = ESMF_regrid(u, Opt) ; Do the regridding for u ;---Use the just generated weight file v_reclin = ESMF_regrid_with_weights(v,Opt@WgtFileName,False) ;---Perform the regrid for wr: WRF ==> rectilinear (_reclin) wr_reclin = ESMF_regrid(wr, Opt) ;---Print rectilinear variable information printVarSummary(u_reclin) print("u_reclin: min="+min(u_reclin)+" max="+max(u_reclin)) print("----") printVarSummary(v_reclin) print("v_reclin: min="+min(v_reclin)+" max="+max(v_reclin)) print("----") printVarSummary(wr_reclin) print("wr_reclin: min="+min(wr_reclin)+" max="+max(wr_reclin)) print("----") nmsg = num(ismissing(u_reclin)) print("RectilinearGrid: nmsg="+nmsg) print("----") ;exit() ;====================================================================== ;Use 'advect_variable' on the rectilinear grid ;====================================================================== adv_vrt = advect_variable(u_reclin, v_reclin, wr_reclin, 0, "vorticity advection", "(m/s)(s-1)",0) printVarSummary(adv_vrt) printMinMax(adv_vrt,0) print("----") ;exit() ;====================================================================== ; interpolate the advected quantity on the rectilinear grid back onto the original WRF grid ;====================================================================== ;---For clarity, delete the above options and start again. delete(Opt) ;---Options for regridding rectilinear to WRF (curvilinear) grid Opt = True ; The grid descriptions have been generated in 'Part A' ; but we still need to provide their names. Opt@SkipSrcGrid = True Opt@SkipDstGrid = True Opt@DstFileName = "WRF.SCRIP_grid_description."+nlat+"_x_"+nlon+".nc" ; Name of source and Opt@SrcFileName = "Rectilinear.SCRIP_grid_description.nc" ; destination files ; Opt@SrcFileName = "WRF.SCRIP_grid_description."+nlat+"_x_"+nlon+".nc" ; destination files Opt@ForceOverwrite = True Opt@SrcTitle = Opt@SrcFileName ; source grid Opt@SrcMask2D = where(ismissing(u_reclin),0,1) Opt@SrcRegional = True Opt@DstGridType = "curvilinear" Opt@DstTitle = "Rectilinear_to_WRF" Opt@DstRegional = True Opt@InterpMethod = InterpMethod Opt@WgtFileName = "Rect_to_WRF.WgtFile_"+InterpMethod+"."+nlat+"_x_"+nlon+".nc" adv_wrf = ESMF_regrid(adv_vrt,Opt) adv_wrf@long_name = "Advection Vorticity: WRF Grid" printVarSummary(adv_wrf) print("adv_wrf: min="+min(adv_wrf)+" max="+max(adv_wrf)) ;================================================================= ;--- Plotting ;================================================================= level = 50 ; arbitrary ... pick a level ;wks = gsn_open_wks("png","d04_div_"+level+"_"+x) wks = gsn_open_wks("png","/home/lara/Documents/Simulaciones_Ubuntu/STC16/NCL_plots/VorticityAdvection/ESMF_regrid") ;================================================================= ;--- Plot: contour wit map using WRF information ;================================================================= res = True ;;res@gsnMaximize = True ; maximize plot in frame res@gsnDraw = False ; will panel later res@gsnFrame = False res@tiMainFont = "helvetica" ; default is helvetica-bold res@gsnAddCyclic = False res@cnFillOn = True res@cnFillMode = "RasterFill" res@cnFillPalette = "ViBlGrWhYeOrRe" res@cnLinesOn = False res@cnLineLabelsOn= False symMinMaxPlt (adv_wrf(level,:,:),20,False,res) res@cnLevelSpacingF = 0.5*res@cnLevelSpacingF ; for 'fun' ... modify ; res@lbOrientation = "Vertical" res@lbLabelBarOn = False ; will add in panel res@pmTickMarkDisplayMode = "Always" ; nicer tickmarks res@mpDataBaseVersion = "MediumRes" ; better and more map outlines res@mpFillOn = False ;---------------------------------------------------------------------- ; Set resources to zoom in on lat/lon area of interest using either ; the WRF map projection defined on the WRF file, or by setting the ; lat/lon limits ourselves. ;---------------------------------------------------------------------- USE_WRF_PROJECTION = True if(USE_WRF_PROJECTION) then res = wrf_map_resources(a,res) ;---Change some of the defaults set by wrf_map_resources. res@mpOutlineBoundarySets = "AllBoundaries" ; more outlines res@mpDataSetName = "Earth..2" res@mpGeophysicalLineColor = "black" ; gray res@mpNationalLineColor = "Black" ; gray res@mpNationalLineThicknessF = 2.0 ; 0.5 res@mpGeophysicalLineThicknessF = 2.0 ; 0.5 else res@mpProjection = "CylindricalEquidistant" res@mpMinLatF = minLat2D res@mpMaxLatF = maxLat2D res@mpMinLonF = minLon2D res@mpMaxLonF = maxLon2D res@tiMainOffsetYF = -0.03 ; Move the title down end if plot_regrid = gsn_csm_contour_map(wks,adv_wrf(level,:,:),res) end