; 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 from point A to point B ; Add some label info to the Y-axis load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" begin ; ;---Read wrf file ;a = addfile("../../wrfout_d01_2015-09-08_12:00:00","r") ; small domain ouput 5 a = addfile("../wrfout_d01_2015-09-08_12:00:00_big2","r") ; big domain ouput 10 ;a = addfile("./wrfout_d04_2015-09-09_00:00:00","r") ; 1.11km ouput 10 ; read need data xlat = wrf_user_getvar(a, "XLAT",0) xlon = wrf_user_getvar(a, "XLONG",0) ter = wrf_user_getvar(a, "HGT",0) ;--Set up workstation type = "png" ; type = "pdf" ; type = "ps" ; type = "ncgm" type@wkWidth = 3000 ; Increase size for a slightly type@wkHeight = 3000 ; Output folder gridsize="1.11" direction = "WE" folder="./clean_vertical_wind_" + gridsize + "km_" + direction + "_longWE" mkdir = systemfunc ("mkdir -p " + folder) wks = gsn_open_wks(type,folder+"/plt_CrossSection_3") ; set bounding box ; A'->A test ---WE--- ;lats = (/ 36.59845, 36.585503 /) ;lons = (/ 139.283992, 138.602908 /) ;xoption=0 ; B->B' test ---SN--- ;lats = (/ 34.946092, 36.7485 /) ;lons = (/ 139.028067, 138.796036 /) ;xoption=1 ;---Please select your plot case if (direction .eq. "WE") then ; C->C' ---WE--- lats = (/ 36.59845, 36.59845 /) lons = (/ 136.50, 142.0 /);(/ 138.50, 141.0 /) ; this is origin xoption=0 else if (gridsize .eq. "30") then ;D->D' ---SN for 30 km--- lats = (/ 34.548864,36.881694 /) lons = (/ 140.13045,139.665556 /) xoption=1 end if if (gridsize .eq. "1.11") then ;E->E' ---SN for 1.11 km--- lats = (/ 34.902453,36.751536 /) lons = (/ 139.090669,138.798081 /) xoption=1 end if end if ; X-axis label ;xoption=0 ;0 means lon X-axis, 1 means lat X-axis if (xoption .eq. 0) then X_desc = "longitude" xvalue=xlon px = 0.5;0.2 end if if (xoption .eq. 1) then X_desc = "latitude" xvalue = xlat px = 0.1 end if ;---Set zoomin indices ; loc(0,;) is west-east (x) ; loc(1,:) is south-north (y) ; subtract one since we want to use it as an index in NCL loc = wrf_user_ll_to_ij(a, lons, lats, True) loc = loc - 1 x_start = loc(0,0);lon x_end = loc(0,1) y_start = loc(1,0) y_end = loc(1,1) ;---Set some basic resources res = True res@gsnMaximize =True res@gsnDraw = False res@gsnFrame = False ;---Read time 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) ;-------------------------------------------------------------------------- ; Loop of time to do plot ;--------------------------------------------------------------- ;do it = 356,356;ntimes-1,18 ; for output5 when need 0909_1740 (UTC) 30km do it = 178,178;ntimes-1,18 ; for output10 when need 0909_1740 (UTC) 30km ;do it = 106,106;ntimes-1,18 ; for output10 when need 0909_1740 (UTC) 1.11km print("Working on time: " + times(it) ) rh = wrf_user_getvar(a,"rh",it) ; relative humidity z = wrf_user_getvar(a, "z",it) ; grid point height u = wrf_user_getvar(a,"ua",it) ; u in m/s w = wrf_user_getvar(a,"wa",it) ; w in m/s v = wrf_user_getvar(a,"va",it) ; w in m/s if ( FirstTime ) then ; get height info for labels zmin = 0. zmax = 12;max(z)/1000. nz = floattoint(zmax + 1) FirstTime = False end if ;--------------------------------------------------------------- ; Plot a cross session that run from point A to point B plane = new(4,float) plane = (/ x_start,x_end, y_start,y_end /) ; start x;y & end x;y point opts = True ; start and end points specified rh_plane = wrf_user_intrp3d(rh,z,"v",plane,0.,opts) u_plane = wrf_user_intrp3d(u,z,"v",plane,0.,opts) w_plane = wrf_user_intrp3d(w,z,"v",plane,0.,opts) v_plane = wrf_user_intrp3d(v,z,"v",plane,0.,opts) X_plane = wrf_user_intrp2d(xvalue,plane,0,opts) dim = dimsizes(rh_plane) ; Find the data span - for use in labels zspan = dim(0) ; Options for XY Plots ; X-axis lables dimsX = dimsizes(X_plane) xmin = X_plane(0) xmax = X_plane(dimsX(0)-1) xspan = dimsX(0)-1 nx = floattoint( (xmax-xmin)/px + 1) print ("dimsX= "+dimsX+" xmin= "+xmin+" xmax= "+xmax+" xspan= "+xspan+ " nx="+nx) ;--------------------------------------------------------------- ; Options for XY Plots opts_xy = res opts_xy@tiXAxisString = X_desc opts_xy@tiYAxisString = "Height (km)" opts_xy@tmXTOn = False opts_xy@tmYROn = False opts_xy@tmXBMode = "Explicit" if (nx .lt. 0) then nx = -nx end if opts_xy@tmXBValues = fspan(0,xspan,nx) ; Create tick marks opts_xy@tmXBLabels = sprintf("%.1f",fspan(xmin,xmax,nx)) ; Create labels opts_xy@tmXBLabelFontHeightF = 0.015 opts_xy@tmYLMode = "Explicit" opts_xy@tmYLValues = fspan(0,zspan,nz) ; Create tick marks opts_xy@tmYLLabels = sprintf("%.1f",fspan(zmin,zmax,nz)) ; 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 ; Plotting options for RH opts_rh = opts_xy opts_rh@cnLinesOn = False opts_rh@cnFillOn = True opts_rh@lbOrientation = "vertical" opts_rh@ContourParameters = (/ 10., 90., 10. /) opts_rh@cnFillColors = (/"White","White","White", \ "White","Chartreuse","Green", \ "Green3","Green4", \ "ForestGreen","PaleGreen4"/) mainstring = times(it) + " (UTC)" leftstring = "Relative humidity %" rightstring = "Wind (m/s)" size = 0.02 opts_rh@tiMainString = mainstring opts_rh@gsnLeftString = leftstring ; add the gsn titles opts_rh@gsnRightString = rightstring opts_rh@tiMainFontHeightF = 0.030 opts_rh@gsnLeftStringFontHeightF = size ; instead of using txFontHeightF or gsnStringFontHeightF opts_rh@gsnRightStringFontHeightF = size ; individually set each string's font height. ; Get the contour info for the rh and temp contour_rh = gsn_csm_contour(wks,rh_plane,opts_rh) ;------------------------------------------------ ; vector plot ;------------------------------------------------ vecres = opts_xy ; vector only resources vecres@vcGlyphStyle = "LineArrow" vecres@vcLineArrowThicknessF = 6 vecres@vcMinDistanceF = 0.008 vecres@vcRefLengthF = 0.04 ;0 means lon X-axis: use u, w, 1 means lat X-axis: use v, w if (xoption .eq. 0) then vecres@vcRefMagnitudeF = 10 end if if (xoption .eq. 1) then vecres@vcRefMagnitudeF = 25 end if ; set for overlay vecres@gsnLeftString = " " ; add the gsn titles vecres@gsnRightString = " " ; Plot vector ;0 means lon X-axis: use u, w, 1 means lat X-axis: use v, w if (xoption .eq. 0) then vector = gsn_csm_vector(wks,u_plane,w_plane,vecres) end if if (xoption .eq. 1) then vector = gsn_csm_vector(wks,v_plane,w_plane,vecres) end if ;str_res = True ;str_res@stMinArrowSpacingF = 0.6 ; arrow spacing. ;str_res@stArrowLengthF = 0.1 ; changes the size of the arrows. ;stream_line = gsn_csm_streamline(wks,v_plane,w_plane,vecres) ; MAKE PLOTS overlay(contour_rh,vector) draw(contour_rh) frame(wks) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; end do ; END OF TIME LOOP end