begin
;---Open WRF output file.
idir = "/home/numodel2/SOULIK/WRFOUT/"
ifil = "wrfout_d01_2018-08-15_12:00:00"
a = addfile(idir + ifil, "r")
;---Set some basic plot options
type = "pdf" ; Choose x11/png/pdf
wks = gsn_open_wks(type,"wrf_pressurelevel")
res = True
res@gsnFrame = False
res@gsnDraw = False
res@gsnMaximize = True
res@gsnLeftString = ""
res@gsnRightString = ""
;---Necessary for contours to be overlaid correctly on WRF projection
res@tfDoNDCOverlay = True ; Tell NCL you are doing a native plot
;---The specific pressure levels that we want the data interpolated to.
pressure_levels = (/ 850./) ; pressure levels to plot
nlevels = dimsizes(pressure_levels) ; number of pressure levels
;---Read variable we will need
times = tostring(a->Times)
nt = 300
lat2d = a->XLAT(0,:,:) ; Latitude
lon2d = a->XLONG(0,:,:) ; Longitude
do it = 48, nt-1 ; TIME LOOP
print("Working on time: " + times(it) )
tc = wrf_user_getvar(a,"tc",it) ; T in C
u = wrf_user_getvar(a,"ua",it) ; u averaged to mass points
v = wrf_user_getvar(a,"va",it) ; v averaged to mass points
p = wrf_user_getvar(a, "pressure",it) ; pressure is our vertical coordinate
z = wrf_user_getvar(a, "z",it) ; grid point height
do level = 0,nlevels-1 ; LOOP OVER LEVELS
pressure = pressure_levels(level)
tc_plane = wrf_user_intrp3d(tc,p,"h",pressure,0.,False)
z_plane = wrf_user_intrp3d( z,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)
wrf_smooth_2d( tc_plane, 3 )
wrf_smooth_2d( z_plane, 3 )
; Plotting options for Tc
opts = res
opts@tiMainString = times(it)
opts@gsnLeftString = tc_plane@description+" ("+tc_plane@units+") at "+tc_plane@PlotLevelID+"~C~" + \
z_plane@description+" ("+z_plane@units+") at "+z_plane@PlotLevelID+"~C~" + \
"Wind (" + u_plane@units + ") at "+tc_plane@PlotLevelID
opts@gsnLeftStringFontHeightF = 0.01
opts = wrf_map_resources(a,opts)
opts@mpDataBaseVersion = "MediumRes" ; better and more map outlines
opts@mpDataSetName = "Earth..4"
opts@mpOutlineBoundarySets = "AllBoundaries"
opts@mpOutlineOn = True
opts@mpGeophysicalLineColor = "Black"
opts@cnLineColor = "Red"
opts@cnInfoLabelOn = False
opts@gsnContourLineThicknessesScale = 2.0
contour_tc = gsn_csm_contour_map(wks,tc_plane,opts)
delete(opts)
; Plotting options for Wind Vectors
opts = res
opts@vcMinDistanceF = 0.02
opts@vcRefLengthF = 0.02
opts@vcMinFracLengthF = 0.2
opts@vcGlyphStyle = "WindBarb"
opts@vcRefAnnoOn = False
vector = gsn_csm_vector(wks,u_plane,v_plane,opts)
delete(opts)
; Plotting options for Geopotential Height
opts = res
printMinMax(z_plane,0)
opts@cnLevelSelectionMode = "ManualLevels"
opts@cnMinLevelValF = 1000
if ( pressure .eq. 850 ) then
opts@cnMaxLevelValF = 2000
opts@cnLevelSpacingF = 20
else if ( pressure .eq. 700 ) then
opts@cnMaxLevelValF = 3300
opts@cnLevelSpacingF = 30
else if ( pressure .eq. 500 ) then
opts@cnMaxLevelValF = 5820
opts@cnLevelSpacingF = 60
else if ( pressure .eq. 200 ) then
opts@cnMaxLevelValF = 13200
opts@cnLevelSpacingF = 60
end if
end if
end if
end if
opts@cnLineColor = "Blue"
opts@cnLineLabelBackgroundColor = -1
opts@cnInfoLabelOn = False
opts@gsnContourLineThicknessesScale = 3.0
contour_z = gsn_csm_contour(wks,z_plane,opts)
delete(opts)
;---Overlay plots on map and draw.
overlay(contour_tc,contour_z)
overlay(contour_tc,vector)
draw(contour_tc) ; This will draw all overlaid plots and the map
frame(wks)
end do ; END OF LEVEL LOOP
end do ; END OF TIME LOOP
end
Commentaires