; 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 = 12 ; circle at first gsres@gsMarkerThicknessF =3 gsres@gsMarkerSizeF = 0.021 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 Meteo_Lat = -16.385892 ; SArequipa_Lat = -16.385892;39.8 Meteo_Lon = -71.533015 ; SArequipa_Lon = -71.533015 ; 116.4667 ; 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 = "NARR3dWSM6_CONUS_UCM_YSU_AugMean_noMic" do iday = 14,14; idomain = 2 ;1 files = systemfunc("ls ./wrfout_d0"+idomain+"_20??-??-??_*: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 = "pdf" type@wkWidth = 800 type@wkHeight = 800 plotwinds= True ; for 2021102700 ; plotwinds= False ; if(plotwinds) then figurename = "wrfout_d0"+idomain+"_thetaWinds_crossSN_"+ifiles else figurename = "wrfout_d0"+idomain+"_theta_crossSN_"+ifiles end if wks = gsn_open_wks(type,figurename) gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ; select color map ; gsn_define_colormap(wks,"BlAqGrWh2YeOrReVi22") ; white in middle ; cmap = gsn_retrieve_colormap(wks) ;do icorlor = 2,4 ; cmap(icorlor,:) = cmap(0,:) ;end do ; gsn_define_colormap(wks,cmap) ; 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->so2(it,:,:,:) ; u on mass points temp = (/temp*1000./) SO2total = temp SO2total@units = "ppbv";(/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) SO2total_plane = wrf_user_intrp3d(SO2total,z,"v",plane,angle,opts) dim = dimsizes(theta_plane) ; Find the data span - for use in labels print(dim) resolution = (/15,3.,0.16/) resolution_in_use = resolution(idomain-1) zspan = dim(0) deltava_plane = va_plane deltava_plane(0,:) = (/0./) deltaz = (zmax-zmin)/1000. ; in km deltawa_plane = wa_plane deltawa_plane(:,0) = (/0./) deltawa_plane(:,1:) = (/(wa_plane(:,1:)-wa_plane(:,0:dim(1)-2))/resolution_in_use/) deltava_plane(1:,:) = (/(va_plane(1:,:)-va_plane(0:zspan-2,:))/deltaz/) vorticity = deltava_plane vorticity = (/deltawa_plane-deltava_plane/) print("deltaz="+deltaz+" horizontal resolution_in_use ="+resolution_in_use) 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 = 400 ; domain 1 n_YMin = 100 ; domain 1 ; n_YMax = 50 ; domain 1 ; n_XMax = id_d+5 ; domain 1 n_XMin = 270+5 ; narrower n_XMax = n_XMin + 37-5-7-5 factorxhu = 4 n_XMin = (265+5)*factorxhu ; wider n_XMax = n_XMin + (80-15-10-10-8)*factorxhu; 101 ; 51 opts_xy@trYMaxF = n_YMax opts_xy@trYMinF = n_YMin opts_xy@trXMaxF = n_XMax opts_xy@trXMinF = n_XMin ; 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@cnFillMode = "AreaFill" ; opts_theta@ContourParameters = (/ 302., 314., 1. /) ; opts_theta@ContourParameters = (/305,325,1/); (/ 1.94, 1.97, .0025 /) opts_theta@ContourParameters = (/310,328,1/); for Fig. 6 (/ 1.94, 1.97, .0025 /) ; opts_theta@ContourParameters = (/-10,10,.2/); (/ 1.94, 1.97, .0025 /) 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(-400,400,10)/(resolution(idomain-1)/factorxhu/cos(angle*pi/180)) + id_d opts_theta@tmXBLabels = ispan(-400,400,10) 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 = " SO2~C~ ppb"; +x@units ; "O~B~3, ~N~ppbv" ; opts_theta@lbTitleString = " vorticity,10~S~-3~N~s~S~-1"; +x@units ; "O~B~3, ~N~ppbv" opts_theta@lbTitleString = " ~F8~q~F~, K"; +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 opts_theta@tmXBLabelDeltaF = -0.5 ; opts_theta@tiXAxisOffsetYF = 0.512 ; #does not seem to work ; print(opts_theta) ; contour_theta = wrf_contour(a,wks,WSP_plane,opts_theta) ; contour_SO2 = wrf_contour(a,wks,deltava_plane, opts_theta) ; contour_SO2 = wrf_contour(a,wks,vorticity, opts_theta) contour_SO2 = wrf_contour(a,wks,theta_plane, opts_theta) opts_theta@ContourParameters = (/305,345,1/); (/ 1.94, 2., .005 /) opts_theta@cnFillOn = False ; opts_theta@cnInfoLabelPerimOn = False ; opts_theta@cnInfoLabelOn = False ; opts_theta@cnLineLabelPerimOn = False; opts_theta@cnLineLabelPerimColor = "white" ; opts_theta@cnInfoLabelPerimThicknessF = 0.000001 ; opts_theta@cnLineLabelFontHeightF = 0.0081 ; opts_theta@cnLineColor = "white"; opts_theta@cnLineThicknessF = 4.5 opts_theta@cnLineLabelsOn = True ; contour_theta = wrf_contour(a,wks,theta_plane,opts_theta) opts_theta@cnLineThicknessF = 2.0 opts_theta@cnLineColor = "blue"; opts_theta@cnLineLabelsOn = False ; contour_theta2 = wrf_contour(a,wks,theta_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@trYMinF = n_YMin 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/factorxhu,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@trYMinF = zcordinate(n_YMin)*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+ter_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,opts_ter@trYMinF+85.08,gsres); terrain in m tres = True tres@txFontHeightF = 0.017 tres@txFontColor= "red" gsn_text(wks,contour_ter,"Arequipa", id_d,opts_ter@trYMinF+485.08,tres) ; not working dmmytxt = gsn_add_text(wks,contour_ter,"Arequipa", id_d,opts_ter@trYMinF+585.08,tres) ; working tres@txAngleF = 90 dmmytxt2 = gsn_add_text(wks,contour_ter,"Chachani", id_d+7.5*factorxhu,opts_ter@trYMinF+2685.08-200,tres) ; working tres@txAngleF = 0 Times_char = a->Times tres@txFontColor= "black" if(chartostring(Times_char(0,0:6)).eq."2023-08") then dmmytxtabc = gsn_add_text(wks,contour_ter, "(f)", id_d+20*factorxhu,opts_ter@trYMinF+445.08,tres) ; working else print(chartostring(Times_char(0,0:6))+"not xxxfinding 2023-08") end if if(chartostring(Times_char(0,0:6)).eq."2023-07") then dmmytxtabc = gsn_add_text(wks,contour_ter, "(e)", id_d+20*factorxhu,opts_ter@trYMinF+445.08,tres) ; working end if if(chartostring(Times_char(0,0:6)).eq."2023-05") then dmmytxtabc = gsn_add_text(wks,contour_ter, "(d)", id_d+20*factorxhu,opts_ter@trYMinF+445.08,tres) ; working end if if(chartostring(Times_char(0,0:6)).eq."2022-03") then dmmytxtabc = gsn_add_text(wks,contour_ter, "(c)", id_d+20*factorxhu,opts_ter@trYMinF+445.08,tres) ; working end if if(chartostring(Times_char(0,0:12)).eq."2021-10-28_12") then dmmytxtabc = gsn_add_text(wks,contour_ter, "(b)", id_d+20*factorxhu,opts_ter@trYMinF+445.08,tres) ; working end if if(chartostring(Times_char(0,0:12)).eq."2021-10-28_08") then dmmytxtabc = gsn_add_text(wks,contour_ter, "(a)", id_d+20*factorxhu,opts_ter@trYMinF+445.08,tres) ; working end if if(chartostring(Times_char(0,0:12)).eq."2021-10-28_09") then dmmytxtabc = gsn_add_text(wks,contour_ter, "(b)", id_d+20*factorxhu,opts_ter@trYMinF+445.08,tres) ; working end if if(chartostring(Times_char(0,0:12)).eq."2021-10-28_13") then dmmytxtabc = gsn_add_text(wks,contour_ter, "(c)", id_d+20*factorxhu,opts_ter@trYMinF+445.08,tres) ; working end if if(chartostring(Times_char(0,0:12)).eq."2021-10-28_15") then dmmytxtabc = gsn_add_text(wks,contour_ter, "(d)", id_d+20*factorxhu,opts_ter@trYMinF+445.08,tres) ; working end if if(chartostring(Times_char(0,0:9)).eq."2021-10-19") then dmmytxtabc = gsn_add_text(wks,contour_ter, "(a)", id_d+20*factorxhu,opts_ter@trYMinF+445.08,tres) ; working end if ; 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) LT_hr = mod(stringtoint(chartostring(Times_char(0,11:12)))-5+24, 24) ; plot = wrf_overlays(a,wks,(/contour_SO2,contour_theta,vector,contour_ter/),pltres) if(plotwinds) then plot = wrf_overlays(a,wks,(/contour_SO2,contour_theta,contour_theta2, vector, contour_ter/),pltres) else plot = wrf_overlays(a,wks,(/contour_SO2,contour_theta,contour_theta2, contour_ter/),pltres) end if setvalues contour_SO2 ; "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,:12) + "UTC "+ LT_hr+"LT " "tiMainOffsetYF" : -0.0321 ; "gsnRightString" : "RString" ; res@TimeLabel ; "gsnPaperOrientation" : "portrait" "tmXBLabelsOn" : True "tiXAxisString" : "Distance, km" "tiXAxisOffsetYF" : 0.12 ; work "tmXBLabelDeltaF": -0.5 end setvalues draw(contour_SO2) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; end do ; END OF TIME LOOP delete(ter_plane) delete(wa_plane) delete(tc_plane) delete(theta_plane) delete(va_plane) delete(deltava_plane) delete(vorticity) delete(deltawa_plane) delete(ua_plane) delete(WSP_plane); = wrf_user_intrp3d(WSP,z,"v",plane,angle,opts) delete(zcordinate) delete(SO2total_plane) frame(wks) ; 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("rsync -av "+figurename+".png xhu@cumulus.caps.ou.edu:/var/www/html/micronet/Brazil/Simulations/4.3.3/wrfchem4.3.3FNL025_2d_Peru_WACCMICBC_Sabancaya.2023080600/.") system("mv "+figurename+"."+type+" ./figures/.") end do ; ifiles end do ; iday ; system("rsync -av wrf_CrossSection_SO2plusVecter_16tilt_addTheta_vorticity.ncl xhu@cumulus.caps.ou.edu:/var/www/html/micronet/Brazil/Simulations/4.3.3/wrfchem4.3.3FNL025_2d_Peru_WACCMICBC_Sabancaya.2023080600/.") end