Temperature, Geopotential Height & Wind

---

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