; xhu@ou.edu Feb. 15, 2022
; Example script to produce plots for a WRF real-data run,
; with the ARW coordinate dynamics option.
; Plot data on a cross section
; This script will plot data at a set angle through a specified point
; Add some label info to the Y-axis
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/wrf/WRFUserARW.ncl"
;load "/home1/01178/tg804586/./tools/NCl/6.3.0/WRFUserARW.ncl"
begin
gsres = True
gsres@gsMarkerIndex = 16 ; circle at first
gsres@gsMarkerThicknessF = 1
gsres@gsMarkerSizeF = 0.01
gsres@gsMarkerColor = "black"
gsres@gsMarkerColor = "red"
; gsn_polymarker(wks,plot,-96.8,32.7833333,gsres)
; gsn_text(wks,plot,"D",-96.8,32.7833333+0.1,tres)
; Meteo_Lat = 32.7833333
; Meteo_Lon = -96.8
Meteo_Lat = 35.534614
Meteo_Lon = -98.099210 ; Enlink
; simname = "NARR3dWSM6_CONUS_UCM_YSU"
;do iday = 7,7; 29,30 ; 3,3; 4,10
; simname = "NARR3dWSM6_CONUS_UCM_YSU_JulAugMean"
; simname = "NARR3dWSM6_CONUS_UCM_YSU_JulAugMean_noMic"
simname = "chem4.3.3LES3d_Hu2021JGR_CH4NEI2017_Wetchart131_agwasteOce"
; simname = "NARR3dWSM6_CONUS_UCM_YSU_AugMean_noMic"
do iday = 14,14;
if (iday.lt.10) then
folder = "wrf"+simname+".2022020"+iday+"00"
else
folder = "wrf"+simname+".202202"+iday+"18"
end if
idomain = 3 ;1
files = systemfunc("ls file_target_to_be_substitute/wrfout_d0"+idomain+"_202?-??-??_*:00:00")
do ifiles = 0 , dimsizes(files)-1
;do ifiles = 1,1; 0 , dimsizes(files)-1
a = addfile(files(ifiles)+".nc","r")
opt = True
ij = wrf_user_ll_to_ij(a,Meteo_Lon, Meteo_Lat,opt)
print("ij = "+ij(0)+","+ij(1))
; We generate plots, but what kind do we prefer?
type = "png"
type@wkWidth = 800
type@wkHeight = 800
figurename = "wrfout_d0"+idomain+"_CH4_crossSN_"+ifiles
wks = gsn_open_wks(type,figurename)
gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ; select color map
; Set some basic resources
res = True
; res@MainTitle = "REAL-TIME WRF"
res@Footer = False
ter_res = True
opts_ter = ter_res
opts_ter@gsnDraw = False
opts_ter@gsnFrame = False
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
times = wrf_user_getvar(a,"times",-1) ; get times in the file
ntimes = dimsizes(times) ; number of times in the file
FirstTime = True
mdims = getfilevardimsizes(a,"P") ; get some dimension sizes for the file
nd = dimsizes(mdims)
;---------------------------------------------------------------
do it = 0,ntimes-1,2 ; TIME LOOP
print("Working on time: " + times(it) )
res@TimeLabel = times(it) ; Set Valid time to use on plots
ter = wrf_user_getvar(a,"ter",0)
xlon = wrf_user_getvar(a, "XLONG",0)
xlat = wrf_user_getvar(a, "XLAT",0)
tc = wrf_user_getvar(a,"tc",it) ; T in C
theta = wrf_user_getvar(a,"theta",it) ; relative humidity
z = wrf_user_getvar(a, "z",it) ; grid point height
va = wrf_user_getvar(a,"va",it) ; u on mass points
ua = wrf_user_getvar(a,"ua",it) ; u on mass points
WSP = sqrt(ua^2+va^2)
wa = wrf_user_getvar(a,"wa",it) ; v on mass points
PBLH = a->PBLH(0,:,:)
temp = a->CH4_TST(it,:,:,:) ; u on mass points
temp = (/a->CH4_TST(0,:,:,:) + a->CH4_ANT(0,:,:,:) - a->CH4_BCK(0,:,:,:)/)
CH4total = temp
CH4total@units = "ppmv";(/temp*1000/)
if ( FirstTime ) then ; get height info for labels
zmin = 0.
zmax = max(z)/1000.
nz = floattoint(zmax/2 + 1)
FirstTime = False
; print("will use "+nz +" layers")
end if
z_values = fspan(zmin,zmax,nz)
;printVarSummary(z_values)
dimzz = dimsizes(z_values)
;---------------------------------------------------------------
; Plot a cross session that run south-north through the middle of the plot
; For this we need a pivot point and a angle
; |
; angle=0 is |
; angle = 90 ; E-W
; angle = -16 ;
angle = 0 ; S-N
plane = new(2,float)
plane = (/ ij(0), ij(1) /) ; pivot point is center of domain (x,y)
opts = False ; start and end points not specified
theta_plane = wrf_user_intrp3d(theta,z,"v",plane,angle,opts)
tc_plane = wrf_user_intrp3d(tc,z,"v",plane,angle,opts)
WSP_plane = wrf_user_intrp3d(WSP,z,"v",plane,angle,opts)
va_plane = wrf_user_intrp3d(va,z,"v",plane,angle,opts)
ua_plane = wrf_user_intrp3d(ua,z,"v",plane,angle,opts)
wa_plane = wrf_user_intrp3d(wa,z,"v",plane,angle,opts)
ter_plane = wrf_user_intrp2d(ter,plane,angle,opts)
lat_plane = wrf_user_intrp2d(xlat,plane,angle,opts)
lon_plane = wrf_user_intrp2d(xlon,plane,angle,opts)
PBLH_plane = wrf_user_intrp2d(PBLH,plane,angle,opts)
CH4total_plane = wrf_user_intrp3d(CH4total,z,"v",plane,angle,opts)
dim = dimsizes(theta_plane) ; Find the data span - for use in labels
zspan = dim(0)
zcordinate = fspan(zmin,zmax,zspan)
distance = (lon_plane( :)- Meteo_Lon)^2+(lat_plane(:)-Meteo_Lat)^2
id_d = minind(distance)
print("found grid point lcose to Meteo_Lon/Meteo_Lat at "+id_d)
dim = dimsizes(theta_plane) ; Find the data span - for use in labels
zspan = dim(0)
; print(zspan)
; Options for XY Plots
opts_xy = res
opts_xy@tiYAxisString = "Height (km)"
opts_xy@cnMissingValPerimOn = False ; True
opts_xy@cnMissingValFillColor = 0
opts_xy@cnMissingValFillPattern = 11
opts_xy@tmYLMode = "Explicit"
opts_xy@tmYLValues = fspan(0,zspan,nz*2*2) ; Create tick marks
opts_xy@tmYLLabels = sprintf("%.1f",fspan(zmin,zmax,nz*2*2)) ; Create labels
opts_xy@tiXAxisFontHeightF = 0.020
opts_xy@tiYAxisFontHeightF = 0.020
opts_xy@tmXBMajorLengthF = 0.02
opts_xy@tmYLMajorLengthF = 0.02
opts_xy@tmYLLabelFontHeightF = 0.015
opts_xy@PlotOrientation = tc_plane@Orientation
; n_YMax = 150 ; domain 1
n_YMax = 50 ; domain 1
; n_XMax = id_d+5 ; domain 1
; n_XMin = 15 ; domain 1
n_XMax = 149; id_d+1 ; domain2
n_XMin = 0 ; 25 ; domain 2
factorxhu = 4 ; matching the number in NCL libraray
opts_xy@trYMaxF = n_YMax
opts_xy@trXMaxF = n_XMax*factorxhu
opts_xy@trXMinF = n_XMin*factorxhu
; opts_xy@gsnOrientation = "portrait"
opts_xy@gsnDraw = True ; Forces the plot to be drawn
opts_xy@gsnFrame = True ; Frame advance
opts_xy@tiMainString = "" ;chartostring(a->Times)
opts_xy@vpHeightF = 0.35
opts_xy@vpWidthF = 0.5
; opts_xy@gsnRightString = chartostring(a->Times)
; Plotting options for RH
opts_theta = opts_xy
opts_theta@cnFillMode = "RasterFill"
; opts_theta@ContourParameters = (/ 302., 314., 1. /)
if(ifiles.eq.1) then
opts_theta@ContourParameters = (/ 1.94, 1.97, .0025 /)
else
opts_theta@ContourParameters = (/ 1.94, 2., .005 /)
end if
; opts_theta@ContourParameters = 5.
opts_theta@tmXBFormat = "f"
opts_theta@pmLabelBarOrthogonalPosF = -0.07
opts_theta@cnFillOn = True
opts_theta@gsnSpreadColors = True ; use full range of colormap
; Plotting options for Temperature
opts_tc = opts_xy
opts_tc@cnInfoLabelOrthogonalPosF = 0.00
opts_tc@ContourParameters = (/ 5. /)
; resolution = (/12,4/)
resolution = (/12,0.8,0.16/)
opts_theta@tmXBMode = "Explicit"
pi = 3.1415926535
opts_theta@tmXBValues = ispan(-50,90,5)/(resolution(idomain-1)/cos(angle*pi/180)) + id_d
opts_theta@tmXBLabels = ispan(-50,90,5)
opts_theta@lbTopMarginF = 0.00516
opts_theta@lbBottomMarginF = 0.35; 54
opts_theta@lbTitleFontHeightF = 0.02075 ; make title smaller
opts_theta@lbLabelFontHeightF = 0.01875
opts_theta@tmXBMajorLengthF = 0.02
opts_theta@tmYLMajorLengthF = 0.02
opts_theta@lbTitleOn = True ; turn on title
opts_theta@lbTitlePosition = "Right" ; title position
opts_theta@lbTitleDirection = "Across" ; title direction
opts_theta@lbTitleOffsetF = -0.034
opts_theta@lbTitleJust = "CenterCenter"
opts_theta@lbTitleString = " CH4~C~ ppm"; +x@units ; "O~B~3, ~N~ppbv"
; Get the contour info for the theta and temp
; printVarSummary(tc_plane)
; contour_tc = wrf_contour(a,wks,tc_plane,opts_tc)
opts_theta@gsnMaximize = False
; print(opts_theta)
; contour_theta = wrf_contour(a,wks,WSP_plane,opts_theta)
; contour_theta = wrf_contour(a,wks,theta_plane,opts_theta)
contour_theta = wrf_contour(a,wks,CH4total_plane,opts_theta)
resvector = True
resvector@vpHeightF = 0.35
resvector@vpWidthF = 0.5
resvector@gsnDraw = True ; Forces the plot to be drawn
resvector@gsnFrame = True ; Frame advance
resvector@vcGlyphStyle = "CurlyVector"
resvector@tmXBFormat = "f"
; resvector@trXMinF = 65.
resvector@trYMaxF = n_YMax
resvector@trXMaxF = n_XMax
resvector@trXMinF = n_XMin
resvector@vcRefAnnoOn = True
resvector@vcRefMagnitudeF = 8. ; define vector ref mag
resvector@vcRefLengthF = 0.025 ; define length of vec ref
resvector@vcRefAnnoOrthogonalPosF = -1.0751 ; move ref vector
resvector@vcRefAnnoString2On = False
resvector@vcRefAnnoParallelPosF = 0.9775 ; move ref vector
resvector@vcMinDistanceF = 0.025 ; larger means sparser
resvector@vcLineArrowHeadMaxSizeF = 0.0075 ; default: 0.05 (LineArrow), 0.012 (CurlyVector)
resvector@gsnMaximize = False
; vector = wrf_vector(a,wks,va_plane,wa_plane*100,resvector)
; vector = wrf_vector(a,wks,va_plane*cos(13*pi/180)-ua_plane*sin(13*pi/180),wa_plane*100,resvector)
vector = wrf_vector(a,wks,va_plane*cos(13*pi/180)-ua_plane*sin(13*pi/180),wa_plane*10,resvector)
; dummy = gsn_add_polymarker(wks,vector,ij(1),0.08,gsres)
;Contour terrain cross section
; print(opts_ter)
; opts_ter@trXMaxF = lat_plane(n_XMax)
; opts_ter@trXMinF = lat_plane(n_XMin)
opts_ter@trXMaxF = n_XMax
opts_ter@trXMinF = n_XMin
opts_ter@trYMaxF = zcordinate(n_YMax)*1000
opts_ter@vpHeightF = 0.35
opts_ter@vpWidthF = 0.5
opts_ter@xyLineThicknesses = (/8., 5./)
; opts_ter@gsnYRefLine = 0.0
; opts_ter@gsnAboveYRefLineColor = "black"
; contour_ter = gsn_csm_xy(wks,lat_plane,ter_plane,opts_ter)
ter_plane_mean_plusref = new ((/2,dimsizes(ter_plane)/),float)
ter_plane_mean_plusref(0,:) = (/ter_plane/)
; ter_plane_mean_plusref(1,:) = (/PBLH_plane*0+2400/)
ter_plane_mean_plusref(1,:) = (/PBLH_plane/)
print("max PBLH "+max(PBLH_plane))
; contour_ter = gsn_csm_xy(wks,lat_plane,ter_plane_mean_plusref,opts_ter)
contour_ter = gsn_csm_y(wks,ter_plane_mean_plusref,opts_ter)
; dummy1 = gsn_add_polymarker(wks,contour_ter,lat_plane(id_d),0.08,gsres)
dummy1 = gsn_add_polymarker(wks,contour_ter,id_d,0.08,gsres)
; MAKE PLOTS
pltres = True
pltres@NoTitles = True
pltres@CommonTitles = True; False
pltres@PanelPlot = True
pltres@gsnPaperOrientation = "portrait"
pltres@FramePlot = False
; plot = wrf_overlays(a,wks,(/contour_theta,vector/),pltres)
Times_char = a->Times
CST_hr = mod(stringtoint(chartostring(Times_char(0,11:12)))-6+24, 24)
plot = wrf_overlays(a,wks,(/contour_theta,vector,contour_ter/),pltres)
setvalues contour_theta
; "vpHeightF" : 0.35
; "vpYF" : 0.9-ifile_loop*0.39
; "vpWidthF" : 0.35
; "vpXF" : 0.15+ispecies*0.355
; "tiMainString" : res@TimeLabel + " UTC";"Cross section"
"tiMainString" : Times_char(0,6:12) + "UTC "+ CST_hr+"CST "
"tiMainOffsetYF" : -0.01
; "gsnRightString" : "RString" ; res@TimeLabel
; "gsnPaperOrientation" : "portrait"
"tmXBLabelsOn" : True
"tiXAxisString" : "Distance, km"
end setvalues
draw(contour_theta)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
end do ; END OF TIME LOOP
delete(ter_plane)
delete(wa_plane)
delete(tc_plane)
delete(theta_plane)
delete(va_plane)
delete(ua_plane)
delete(WSP_plane); = wrf_user_intrp3d(WSP,z,"v",plane,angle,opts)
delete(zcordinate)
delete(CH4total_plane)
frame(wks)
system("mv "+figurename+".png outputdirectory/.")
; system("convert -density 300 -rotate -90 -resize 1000x1000 -trim "+figurename+".eps /nsftor/xhu/public_html/HurricaneImpactonO3/WRFV3.6.1/YSU/"+folder+"/"+figurename+".png")
; system("convert -density 300 -resize 1000x1000 -trim "+figurename+".eps /nsftor/xhu/public_html/HurricaneImpactonO3/WRFV3.6.1/YSU/wrfNARR3dWSM6_CONUS_UCM_YSU.2022020700/"+figurename+".png")
system("rm "+figurename+".eps")
end do ; ifiles
print("finish "+folder)
end do ; iday
end