; http://www2.mmm.ucar.edu/wrf/OnLineTutorial/Graphics/NCL/Examples/CROSS_SECTION/wrf_CrossSection_smooth4.htm --> ; Ref: http://www2.mmm.ucar.edu/wrf/OnLineTutorial/Graphics/NCL/Examples/CROSS_SECTION/wrf_CrossSection_add_smooth_terrain4.ncl ; 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 ; This script adds lon/lat info along X-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" ;---------------------------------------------------------------------- undef("add_outlines_to_map") procedure add_outlines_to_map(wks,plot) local fnames, colors, lnres, n begin ; ; These two shapefiles contain waterways and administrative ; geographic info. We will add them to the existing map using ; polylines. ; fnames = (/"./shp_ostan/iran_province"/) + ".shp" colors = (/"Brown"/) ; water, administrative lnres = True ; resources for polylines lnres@gsLineThicknessF = 2.0 ; 2x as thick ; Ehsan: 3.0 ;***Following loop is unnecessary do n=0,dimsizes(fnames)-1 dumstr = unique_string("primitive") lnres@gsLineColor = colors(n) plot@$dumstr$ = gsn_add_shapefile_polylines(wks, plot, fnames(n), lnres) end do end ;---------------------------------------------------------------------- ; Add markers for cities to given map plot. ;---------------------------------------------------------------------- undef("add_places_to_map") procedure add_places_to_map(wks,plot) local f, stids, stnames, sticaos, stlats, stlons, stelvs, stdims, xres, mkres begin ;---Read shapefile with "places" information. f = addfile("./shp_stations/stations1.shp","r") stids = f->stid stnames = f->stnm sticaos = f->ICAOC stlons = f->stlon stlats = f->stlat stelvs = f->stelv stdims = dimsizes(stnames) ;---Set up text resources to label random areas of interest. txres = True txres@txFontHeightF = 0.010 txres@txJust = "topcenter" ; "topright" ;---Set up marker resources to mark random areas of interest. mkres = True mkres@gsMarkerIndex = 16 mkres@gsMarkerColor = "brown" mkres@gsMarkerSizeF = 7. icao_of_interest = (/ "OITR","OITT","OITL","OICS","OITZ","OIGG", \ "OICC","OIHH","OIIK","OICI","OICK","MAAA", \ "ALKK","MASA","OIII","OIAW","OIFS","OIFM", \ "OIIS","OING","OISY","OIBB","OISS","OIYY", \ "OIMB","OIMM","OIMN","OIKB","OIKK","OIZH" /) ;---Looking for interested icao codes icao_of_interest@_FillValue = "missing" do i=0,stdims-1 if(any(sticaos(i).eq.icao_of_interest)) then ii = ind(icao_of_interest.eq.sticaos(i)) icao_of_interest(ii) = "missing" dumstr = unique_string("primitive") ;---Attach the marker to the map. plot@$dumstr$ = gsn_add_polymarker(wks, plot, stlons(i), stlats(i), mkres) dumstr = unique_string("primitive") ; ;---Attach the text string to the map. plot@$dumstr$ = gsn_add_text(wks, plot, " " + sticaos(i), \ stlons(i), stlats(i), txres) delete(ii) end if end do end ;---------------------------------------------------------------------- ;main code begin ; ; The WRF ARW input file. ; This needs to have a ".nc" appended, so just do it. ; a = addfile("../wrfout_d01_2000-01-24_12:00:00.nc","r") ; a = addfile("./wrfout_d02_2018-08-17_12:00:00"+".nc","r") a = addfile("./wrfout_d01_2018-08-22_00:00:00"+".nc","r") ; We generate plots, but what kind do we prefer? ; type = "x11" type = "pdf" ; type = "ps" ; type = "ncgm" wks = gsn_open_wks(type,"plt_CrossSection_smooth4") ; Set some basic resources res = True res@MainTitle = "REAL-TIME WRF" res@Footer = False pltres = True ter_res = True opts_ter = ter_res opts_ter@gsnYRefLine = 0.0 opts_ter@gsnAboveYRefLineColor = "black" opts_ter@gsnDraw = False opts_ter@gsnFrame = False ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FirstTime = True FirstTimeMap = True times = wrf_user_getvar(a,"times",-1) ; get times in the file ntimes = dimsizes(times) ; number of times in the file mdims = getfilevardimsizes(a,"P") ; get some dimension sizes for the file nd = dimsizes(mdims) xlat = wrf_user_getvar(a, "XLAT",0) xlon = wrf_user_getvar(a, "XLONG",0) ter = wrf_user_getvar(a, "HGT",0) ;--------------------------------------------------------------- lat1 = 35.69 ; Tehran (Mehrabad Airport) lon1 = 51.31 lat2 = 34.35 ; Kermanshah lon2 = 47.15 ;--------------------------------------------------------------- optll = True loc = wrf_user_ll_to_ij(a,(/lon1,lon2/),(/lat1,lat2/),optll) latlon = wrf_user_ij_to_ll(a,loc(0,:),loc(1,:),optll) print(loc) print(latlon) ;--------------------------------------------------------------- ; do it = 0,ntimes-1 ; TIME LOOP do it = 0,18, 4 print("Working on time: " + times(it) ) res@TimeLabel = times(it) ; Set Valid time to use on plots tc = wrf_user_getvar(a,"tc",it) ; T in C rh = wrf_user_getvar(a,"rh",it) ; relative humidity z = wrf_user_getvar(a, "z",it) ; grid point height if ( FirstTime ) then ; get height info for labels zmin = 0. zmax = 18. ; We are only interested in the first 18km nz = floattoint(zmax + 1) end if ;opts_ter@trYMaxF = z_values(dimzz-2)*1000 opts_ter@trYMaxF = zmax*1000 ;--------------------------------------------------------------- ; do ip = 1, 4 ; we are doing 3 plots ; all with the pivot point (plane) in the center of the domain ; at angles 0, 45 and 90 ; ; | ; angle=0 is | ; | ; opts = True ; specifying start and end points plane = new(4,float) ; plane = (/ lon1,lat1 , lon2,lat2 /) ; start x;y & end x;y point plane = (/ loc(0,1),loc(1,1) , loc(0,0),loc(1,0) /) ; start x;y & end x;y point ; plane = new(2,float) ; plane = (/ mdims(nd-1)/2, mdims(nd-2)/2 /) ; pivot point is center of domain (x,y) ; opts = False ; if(ip .eq. 1) then ; angle = 90. ; X_plane = wrf_user_intrp2d(xlon,plane,angle,opts) ; X_desc = "longitude" ; end if ; if(ip .eq. 2) then angle = 0. ; X_plane = wrf_user_intrp2d(xlat,plane,angle,opts) ; X_desc = "latitude" ; end if ; if(ip .eq. 3) then ; angle = 45. X_plane = wrf_user_intrp2d(xlon,plane,angle,opts) X_desc = "longitude" ; end if ; if(ip .eq. 4) then ; angle = 20. ; X_plane = wrf_user_intrp2d(xlon,plane,angle,opts) ; X_desc = "longitude" ; end if rh_plane = wrf_user_intrp3d(rh,z,"v",plane,angle,opts) tc_plane = wrf_user_intrp3d(tc,z,"v",plane,angle,opts) ter_plane = wrf_user_intrp2d(ter,plane,angle,opts) print("Max terrain height in plot " + max(ter_plane)) rh_plane2 = rh_plane tc_plane2 = tc_plane cross_dims = dimsizes(rh_plane2) rank = dimsizes(cross_dims) ;printVarSummary(rh_plane2) iz_do = 25 ;what does it mean? do iz = 0,24 ;what do these loop do? iz_do = iz_do-1 do ix = 0,cross_dims(rank-1)-1 if ( ismissing(rh_plane2(iz_do,ix)) ) then rh_plane2(iz_do,ix) = rh_plane2(iz_do+1,ix) end if if ( ismissing(tc_plane2(iz_do,ix)) ) then tc_plane2(iz_do,ix) = tc_plane2(iz_do+1,ix) end if end do end do ; Find the index where 18km is - only need to do this once if ( FirstTime ) then zz = wrf_user_intrp3d(z,z,"v",plane,angle,opts) b = ind(zz(:,0) .gt. zmax*1000. ) zmax_pos = b(0) - 1 if ( abs(zz(zmax_pos,0)-zmax*1000.) .lt. abs(zz(zmax_pos+1,0)-zmax*1000.) ) then zspan = b(0) - 1 else zspan = b(0) end if delete(zz) delete(b) FirstTime = False print(zspan);;; end if ; 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)/2 + 1) ;--------------------------------------------------------------- ; Options for XY Plots opts_xy = res opts_xy@tiXAxisString = X_desc opts_xy@tiYAxisString = "Height (km)" opts_xy@cnMissingValPerimOn = True opts_xy@cnMissingValFillColor = 0 opts_xy@cnMissingValFillPattern = 11 opts_xy@tmXTOn = False opts_xy@tmYROn = False opts_xy@tmXBMode = "Explicit" 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 opts_xy@PlotOrientation = tc_plane@Orientation ; Plotting options for RH opts_rh = opts_xy opts_rh@ContourParameters = (/ 10., 90., 10. /) ; an array of three elements that represent the minimum level, the maximum level, and the level spacing opts_rh@pmLabelBarOrthogonalPosF = -0.1 opts_rh@cnFillOn = True opts_rh@cnFillColors = (/"White","White","White", \ "White","Chartreuse","Green", \ "Green3","Green4", \ "ForestGreen","PaleGreen4"/) ; Plotting options for Temperature opts_tc = opts_xy opts_tc@cnInfoLabelZone = 1 opts_tc@cnInfoLabelSide = "Top" opts_tc@cnInfoLabelPerimOn = True opts_tc@cnInfoLabelOrthogonalPosF = -0.00005 opts_tc@ContourParameters = (/ 5. /) ; contour level spacing ;Contour terrain cross section contour_ter = gsn_csm_xy(wks,X_plane,ter_plane,opts_ter) ; Get the contour info for the rh and temp contour_tc = wrf_contour(a,wks,tc_plane(0:zmax_pos,:),opts_tc) contour_rh = wrf_contour(a,wks,rh_plane(0:zmax_pos,:),opts_rh) contour_tc2 = wrf_contour(a,wks,tc_plane2(0:zmax_pos,:),opts_tc) contour_rh2 = wrf_contour(a,wks,rh_plane2(0:zmax_pos,:),opts_rh) ;--------------------------------------------------------------- ; MAKE PLOTS if (FirstTimeMap) then lat_plane = wrf_user_intrp2d(xlat,plane,angle,opts) lon_plane = wrf_user_intrp2d(xlon,plane,angle,opts) mpres = True ; mpres@gsnFrame = False mpres@mpOutlineBoundarySets = "National" ;---Boundaries mpres@mpDataSetName = "Earth..4" ; 1 through 4 mpres@mpNationalLineColor = "black" ; "mpCountyLineColor" doesn't work mpres@mpNationalLineThicknessF = 2.5 ; "mpCountyLineThicknessF" doesn't work pltres = True pltres@FramePlot = False optsM = res optsM@NoHeaderFooter = True optsM@cnFillOn = True optsM@lbTitleOn = False optsM@bBoxEndCapStyle = "TriangleBothEnds" ; warning:bBoxEndCapStyle is not a valid resource in plt_CrossSection_smooth4_contour at this time contour = wrf_contour(a,wks,ter,optsM) add_outlines_to_map(wks,contour) add_places_to_map(wks,contour) plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres) lnres = True lnres@gsLineThicknessF = 3.0 lnres@gsLineColor = "Red" do ii = 0,dimsX(0)-2 gsn_polyline(wks,plot,(/lon_plane(ii),lon_plane(ii+1)/),(/lat_plane(ii),lat_plane(ii+1)/),lnres) end do frame(wks) delete(lon_plane) delete(lat_plane) pltres@FramePlot = True end if ;plot = wrf_overlays(a,wks,(/contour_rh,contour_tc/),pltres) ; plot x-section ;plot = wrf_overlays(a,wks,(/contour_rh,contour_tc,contour_ter/),pltres) ; plot x-section plot = wrf_overlays(a,wks,(/contour_rh2,contour_tc2,contour_ter/),pltres) ; plot x-section ; Delete options and fields, so we don't have carry over delete(opts_xy) delete(opts_tc) delete(opts_rh) delete(tc_plane) delete(rh_plane) delete(tc_plane2) delete(rh_plane2) delete(X_plane) delete(ter_plane) ; end do ; make next cross section ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FirstTimeMap = False end do ; END OF TIME LOOP end