; 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/WRF_contributed.ncl" load "/home/xhu2/tools/nclscripts/wrf/WRFUserARW.ncl" begin elevation_array = readAsciiTable("ACTAMERICA-Elevation_B200_20160804_R0.ict", 4, "float", 36) elevation_array@_FillValue = -9999. ; Altitude_AGL_m = elevation_array(:,3) average_points = 30; 20 ; Altitude_GPS = elevation_array(::average_points,1) Altitude_time = elevation_array(::average_points,0) average_points_more = 100 Altitude_GPS_sparse = elevation_array(::average_points_more,1)/1000. Altitude_time_sparse = elevation_array(::average_points_more,0) hours_since_midnight_UTC_simple = round(elevation_array(::average_points,0)/3600.*10, 3)/10. hours_since_midnight_UTC_simple_sparse = round(elevation_array(::average_points_more,0)/3600.*10, 3)/10. hours_since_midnight_UTC_simple_sparse_2d = round(elevation_array(::average_points_more,0)/3600.*100, 3)/100. delete(elevation_array) file_data = "ACTAMERICA-PDS_B200_20160804_R1.ICT" if (file_data.eq."ACTAMERICA-PDS_B200_20160711_R0.ict".or.file_data.eq."ACTAMERICA-PDS_B200_20160715_R1.ict".or.file_data.eq."ACTAMERICA-Hskping_c130_20160721_R0.ict") then data_array = readAsciiTable(file_data, 33, "float", 65) else data_array = readAsciiTable(file_data, 34, "float", 67) end if data_array@_FillValue = -9999. lat = data_array(::average_points,2) lon = data_array(::average_points,3) lat_sparse = data_array(::average_points_more,2) lon_sparse = data_array(::average_points_more,3) ; index_exclude = ind(lon.gt.-75.45) ; index_exclude_sparse = ind((hours_since_midnight_UTC_simple_sparse_2d.lt.18.7.or.hours_since_midnight_UTC_simple_sparse_2d.gt.19.333333).or.Altitude_GPS_sparse.gt.2.3) index_exclude_sparse = ind((hours_since_midnight_UTC_simple_sparse_2d.lt.17.205.or.hours_since_midnight_UTC_simple_sparse_2d.gt.19.37592).or.Altitude_GPS_sparse.gt.2.3) delete(data_array) ; data_array = readAsciiTable("ACTAMERICA-PICARRO_C130_20160725_R0.ict", 6, "float", 38) data_array = readAsciiTable("ACTAMERICA-PICARRO_B200_20160804_R0.ict", 6, "float", 38) data_array@_FillValue = -9999. CO2_obs_raw_temp = data_array(:,3) index_good = ind(.not.ismissing(CO2_obs_raw_temp)) CO2_obs_raw = CO2_obs_raw_temp(index_good) CO2_time_raw = data_array(index_good,2) CO2_obs_raw_runave = runave_Wrap (CO2_obs_raw, average_points, 1) ; 20 point running average with reflective boundary condition ; print(CO2_obs_raw) CO2_obs = linint1 (CO2_time_raw, CO2_obs_raw_runave, False, Altitude_time, 0) ; CO2_obs(index_exclude) = CO2_obs@_FillValue CO2_obs_raw_runave = runave_Wrap (CO2_obs_raw, 120, 1) ; 20 point running average with reflective boundary condition CO2_obs_sparse = linint1 (CO2_time_raw, CO2_obs_raw_runave, False, Altitude_time_sparse, 0) CO2_obs_sparse(index_exclude_sparse) = CO2_obs_sparse@_FillValue dateshow_str = flt2string(hours_since_midnight_UTC_simple_sparse) dateshow_str(index_exclude_sparse) = "" ; dateshow_str@_FillValue print(CO2_time_raw(:10)) print(Altitude_time(:10)) if (dimsizes(Altitude_GPS).ne.dimsizes(CO2_obs)) then print("records in two files doesn't match !!!!!!!!!!!!!!!!" + "data line "+dimsizes(CO2_obs)+" elevation line "+dimsizes(Altitude_GPS)) end if gsres = True gsres@gsMarkerIndex = 16 ; circle at first gsres@gsMarkerThicknessF = 1 gsres@gsMarkerSizeF = 0.01 gsres@gsMarkerColor = "black" gsres@gsMarkerColor = "red" ; 32.453037,-93.834152 Shreveport Regional Airport ; Meteo_Lat = 32.453037 ; Meteo_Lon = -93.834152 ;46.758717,-103.258858 for C130 track on Aug 5 ; Meteo_Lat = 46.758717 ; Meteo_Lon = -103.258858 ; 41.765373,-84.836205 one point along flight starting from NASA ; Meteo_Lat = 41.765373 ; Meteo_Lon = -84.836205 ; 40.856434,-96.763359 Lincoln Airport Nebraska Meteo_Lat = 40.856434 Meteo_Lon = -96.763359 idomain = 1 ; files = systemfunc("ls wrfout_d0"+idomain+"_2011-08-??_*:00:00") ; files = systemfunc("ls wrfout_d01_2016-08-05_1[5-9]* wrfout_d01_2016-08-05_2[0-2]*:00") files = systemfunc("ls wrfout_d01_2016-08-0*00:00 wrfout_d01_2016-08-05_2[0-2]*:00") ; files = systemfunc("ls wrfout_d01_2016-07-25_1[5-9]*00 wrfout_d01_2016-07-25_20*:00") ;do ifiles =18,18; 0 , dimsizes(files)-1 do ifiles = hour_to_be_substitute, hour_to_be_substitute ;do ifiles =19,19; 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) ; We generate plots, but what kind do we prefer? type = "png" type@wkWidth = 1501 type@wkHeight = 1501 ifiles_show = ifiles+ 6 ; figurename = "wrfout_d0"+idomain+"_CO2total_crossSNpartial_overlayObs_combHor_"+ifiles_show ; figurename = "wrfout_d0"+idomain+"_CO2total_crossSNpartial_overlayObs_combHor_"+ifiles ; figurename = "wrfout_d0"+idomain+"_CO2total_crossSNpartial_overlayObs_combHor_wide_"+ifiles ; figurename = "wrfout_d0"+idomain+"_CO2_6components_crossSNpartial_overlayObs_" +ifiles ; figurename = "wrfout_d0"+idomain+"_CO2total_crossSNpartial_overlayObs_combHor_wide_plusT_1"; theta WSP_crossSNpartial = False ; True if (WSP_crossSNpartial) then figurename = "wrfout_d0"+idomain+"_WSP_crossSNpartial_"+ifiles else ; figurename = "wrfout_d0"+idomain+"_T_crossSNpartial_"+ifiles figurename = "wrfout_d0"+idomain+"_QVAPOR_crossSNpartial_"+ifiles end if 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@gsnYRefLine = 0.0 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) PBLH = a->PBLH(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 thetae = wrf_user_getvar(a,"eth",it) ; relative humidity z = wrf_user_getvar(a, "z",it) ; grid point height height = wrf_user_getvar(a,"height",0) 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) EXCH_H = a->EXCH_H(it,1:,:,:) ; u on mass points QVAPOR = a->QVAPOR(it,:,:,:) ; u on mass points QVAPOR = (/QVAPOR*1000./) ; conver to g/kg QCLOUD = a->QCLOUD(it,:,:,:) ; u on mass points QCLOUD = (/a->QCLOUD(it,:,:,:)+a->QGRAUP(0,:,:,:)+a->QICE(0,:,:,:)+a->QRAIN(0,:,:,:)+a->QSNOW(0,:,:,:)/) ; u on mass points QCLOUD = (/QCLOUD*1000.*1000./) ; conver to mg/kg temp = a->CO2_BIO(it,:,:,:) ; u on mass points temp = (/a->CO2_BIO(0,:,:,:) + a->CO2_ANT(0,:,:,:) - a->CO2_BCK(0,:,:,:)/) CO2total = temp CO2total@units = "ppmv";(/temp*1000/) CO2bio = CO2total CO2bio = (/a->CO2_BIO(0,:,:,:) - a->CO2_BCK(0,:,:,:)/) CO2Resp = CO2total CO2Resp = (/a->CO2_BIO_R(0,:,:,:) - a->CO2_BCK(0,:,:,:)/) CO2BCK = CO2total CO2BCK = (/a->CO2_BCK(0,:,:,:)/) CO2Anth = CO2total CO2Anth = (/a->CO2_ANT(0,:,:,:) - a->CO2_BCK(0,:,:,:)/) CO2Photo = CO2total CO2Photo = (/a->CO2_BIO_P(0,:,:,:) - a->CO2_BCK(0,:,:,:)/) wa = wrf_user_getvar(a,"wa",it) ; v on mass points CO2_995wrf = wrf_interp_3d_z(CO2total, height, 995.) u_995wrf = wrf_interp_3d_z(ua,height,995.) v_995wrf = wrf_interp_3d_z(va,height,995.) thetae_995wrf = wrf_interp_3d_z(thetae, height, 995.) theta_995wrf = wrf_interp_3d_z(theta, height, 995.) lat2d = a->XLAT(0,:,:) lon2d = a->XLONG(0,:,:) CO2_995wrf@lat2d = lat2d CO2_995wrf@lon2d = lon2d thetae_995wrf@lat2d = lat2d thetae_995wrf@lon2d = lon2d theta_995wrf@lat2d = lat2d theta_995wrf@lon2d = lon2d u_995wrf@lat2d = lat2d u_995wrf@lon2d = lon2d v_995wrf@lat2d = lat2d v_995wrf@lon2d = lon2d if ( FirstTime ) then ; get height info for labels zmin = 0. zmax = max(z)/1000. nz = floattoint(zmax/2 + 1) FirstTime = False end if ;--------------------------------------------------------------- ; 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 = 90 ; E-W ; angle = -16 ; angle = -8.93 ; angle = 115.; July 25 angle = 119. ; Aug 4 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) thetae_plane = wrf_user_intrp3d(thetae,z,"v",plane,angle,opts) tc_plane = wrf_user_intrp3d(tc,z,"v",plane,angle,opts) CO2total_plane = wrf_user_intrp3d(CO2total,z,"v",plane,angle,opts) CO2bio_plane = wrf_user_intrp3d(CO2bio,z,"v",plane,angle,opts) CO2Anth_plane = wrf_user_intrp3d(CO2Anth,z,"v",plane,angle,opts) CO2BCK_plane = wrf_user_intrp3d(CO2BCK,z,"v",plane,angle,opts) CO2Resp_plane = wrf_user_intrp3d(CO2Resp,z,"v",plane,angle,opts) CO2Photo_plane = wrf_user_intrp3d(CO2Photo,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) WSP_plane = wrf_user_intrp3d(WSP,z,"v",plane,angle,opts) QVAPOR_plane = wrf_user_intrp3d(QVAPOR,z,"v",plane,angle,opts) EXCH_H_plane = wrf_user_intrp3d(EXCH_H,z,"v",plane,angle,opts) QCLOUD_plane = wrf_user_intrp3d(QCLOUD,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) PBLH_plane = wrf_user_intrp2d(PBLH,plane,angle,opts) ; X_plane = wrf_user_intrp2d(xlat,plane,angle,opts) X_plane = wrf_user_intrp2d(xlon,plane,angle,opts) lon_plane = wrf_user_intrp2d(xlon,plane,angle,opts) distance = (lon_plane( :)+96.8)^2+(X_plane(:)-32.7833333)^2 id_d = minind(distance) dim = dimsizes(theta_plane) ; Find the data span - for use in labels print(dim) zspan = dim(0) zcordinate = fspan(zmin,zmax,zspan) ; print(zspan) ; Options for XY Plots opts_xy = res opts_xy@tiYAxisString = "Height, km" ; opts_xy@tiXAxisString = "Latitude, ~S~o~N~N" opts_xy@tiXAxisString = "Longitude, ~S~o~N~W" 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@tmXBMode = "Explicit" lon_plane_round = round(lon_plane*10,3)/10. ; opts_xy@tmXBValues = ispan(0,260,20) ; opts_xy@tmXBLabels = X_plane_round(ispan(0,260,20)) opts_xy@tmXBValues = ispan(0,370,25) opts_xy@tmXBLabels = -lon_plane_round(ispan(0,370,25)) opts_xy@tiXAxisFontHeightF = 0.0108 opts_xy@tiYAxisFontHeightF = 0.0108 opts_xy@tiXAxisOffsetYF = 0.0055 opts_xy@tmXBLabelDeltaF = -0.8 opts_xy@tmXBMajorLengthF = 0.0108 opts_xy@tmYLMajorLengthF = 0.0108 opts_xy@tmYLLabelFontHeightF = 0.0108 opts_xy@tmXBLabelFontHeightF = 0.0108 opts_xy@PlotOrientation = tc_plane@Orientation ; n_YMax = 360 ; domain 1 ; Aug 5 ; n_XMax = 240 ; id_d+1 ; domain 1 ; n_XMin = 130 ; domain 1 ; n_XMin = 0 ; domain 1 n_YMax = 300 ; domain 1 ; July 25 n_XMax = 320; 260 ; id_d+1 ; domain 1 n_XMin = 200 ; domain 1 ; Aug 4 if ( WSP_crossSNpartial ) then ; = True n_YMax = 240; 280 ; domain 1 ; used to run only up to 100 mb else n_YMax = 120; 280 ; domain 1 ; used to run only up to 100 mb end if n_XMax = 270-12; 260 ; id_d+1 ; domain 1 n_XMin = 195-25 ; domain 1 opts_xy@trYMaxF = n_YMax 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.253025 opts_xy@vpWidthF = 0.425 opts_xy@gsnRightString = "Vertical cross-section along the C130 track" ; res@TimeLabel opts_xy@FieldTitle = "Vertical cross-section along the C130 track" ; res@TimeLabel ; Plotting options for CO2total opts_theta = opts_xy opts_theta@cnFillMode = "RasterFill" ; opts_theta@ContourParameters = (/ 390., 406., 1. /) opts_theta@ContourParameters = (/ 396., 408., 1. /) ; opts_theta@ContourParameters = (/ 3., 12., 1. /) ; for WSP ; opts_theta@ContourParameters = (/ 2., 15., 1. /) ; for QVAPOR ; opts_theta@ContourParameters = (/ 295., 306., 1. /) ; for Theta opts_theta@tmXBFormat = "f" opts_theta@pmLabelBarOrthogonalPosF = -0.07 opts_theta@cnFillOn = True opts_theta@gsnSpreadColors = True ; use full range of colormap opts_theta@lbTitleOn = True ; turn on title opts_theta@lbTitleString = " CO~B~2~N~~C~~V10~ppmv" ; "O~B~3, ~N~ppbv" ; title string opts_theta@lbTitlePosition = "Right" ; title position opts_theta@lbTitleFontHeightF= .01348 ; make title smaller opts_theta@lbLabelFontHeightF = 0.01308 opts_theta@lbLabelOffsetF = 0.01 opts_theta@lbTitleDirection = "Across" ; title direction ; opts_theta@lbTopMarginF = -0.051002 ; opts_theta@lbBottomMarginF = 0.21 opts_theta@lbTopMarginF = -0.151002 opts_theta@lbBottomMarginF = 0.321 opts_theta@lbLeftMarginF = 0.1251 opts_theta@lbTitleOffsetF = 0.0004 opts_theta@lbTitleJust = "CenterCenter" ; opts_theta@cnFillColors = (/"White","White","White", \ ; "White","Chartreuse","Green", \ ; "Green3","Green4", \ ; "ForestGreen","PaleGreen4"/) ; Plotting options for Temperature opts_tc = opts_xy opts_tc@cnInfoLabelOrthogonalPosF = 0.00 opts_tc@ContourParameters = (/ 5. /) ; Get the contour info for the theta and temp opts_theta@gsnMaximize = False ; CO2total_plane = (/WSP_plane/) ; CO2total_plane = (/theta_plane/) ; CO2total_plane = (/QVAPOR_plane/) CO2total_plane!0 = "height" CO2total_plane!1 = "lon" CO2total_plane&height = fspan(0,zmax,dim(0) ) CO2total_plane&lon = X_plane va_plane!0 = "height" va_plane!1 = "lon" va_plane&height = fspan(0,zmax,dim(0) ) va_plane&lon = X_plane wa_plane!0 = "height" wa_plane!1 = "lon" wa_plane&height = fspan(0,zmax,dim(0) ) wa_plane&lon = X_plane contour_theta = wrf_contour(a,wks,CO2total_plane,opts_theta) ; contour_theta = wrf_contour(a,wks,WSP_plane,opts_theta) ;;;;;;;;;;;;;;;;;add front symbol;;;;;;;;;;;;;;;;;; gsres = True gsres@gsMarkerThicknessF =6 gsres@gsMarkerSizeF = 0.007081002 gsres@gsMarkerColor = "blue" gsres@gsMarkerIndex = 7 ; Use filled dots for markers. ; dummy_front = gsn_add_polymarker(wks,contour_theta,238.5,6,gsres) txres = True txres@txFontHeightF = 0.012 txres@txBackgroundFillColor = "Transparent" ; text_front = gsn_add_text(wks,contour_theta,"front",238.5,.1 ,txres) ;;;;;;;;;;;;;;;;;add front symbol;;;;;;;;;;;;;;;;;; j_Altitude_GPS = Altitude_GPS/(zmax*1000./dim(0)) x_lon = (lon-min(lon_plane))/((max(lon_plane)-min(lon_plane))/dimsizes(lon_plane)) ; this may cause shift do i_points = 0, dimsizes(x_lon)-1 i_lon = closest_val(lon(i_points),lon_plane) x_lon(i_points) = closest_val(lon(i_points),lon_plane) + (lon(i_points)-lon_plane(i_lon))/(lon_plane(i_lon+1)-lon_plane(i_lon)) end do print("max, min of X_plane"+ min(X_plane)+" "+max(X_plane)+"max, min of lat "+min(lat)+" "+max(lat) ) print("max, min of x_lon"+ min(x_lon)+" "+max(x_lon) ) ; dummy1 = gsn_add_polymarker(wks,contour_theta,x_lon,j_Altitude_GPS,gsres) npts = dimsizes(lat) ; print(contour_theta@contour) ; print(contour_theta) ; getvalues contour_theta@contour ; "cnLevels" : levels ; "cnFillColors" : colors ; end getvalues levels = ispan(396, 408, 1) num_distinct_markers = dimsizes(levels)+1 ; number of distinct markers colors = round(fspan(2,201, num_distinct_markers),3) print(colors) dims=dimsizes(CO2_obs) y_new = new((/num_distinct_markers,dims/),float,-999) x_new = new((/num_distinct_markers,dims/),float,-999) do j = 0, num_distinct_markers-1 if (j.eq.0) then indexes = ind(CO2_obs.lt.levels(0)) end if if (j.eq.num_distinct_markers-1) then indexes = ind(CO2_obs.ge.max(levels)) ; print("color "+j+" indexes="+indexes) end if if (j.gt.0.and.j.lt.num_distinct_markers-1) then indexes = ind(CO2_obs.ge.levels(j-1).and.CO2_obs.lt.levels(j)) end if if (.not.any(ismissing(indexes))) then npts_range = dimsizes(indexes) ; # of points in this range. y_new(j,0:npts_range-1) = j_Altitude_GPS(indexes) x_new(j,0:npts_range-1) = x_lon(indexes) print("color "+j+"npts_range="+npts_range) end if delete(indexes) ; Necessary b/c "indexes" may be a different ; size next time. end do ; j = 0, num_distinct_markers-1 gsres = True gsres@gsMarkerIndex = 4 ; circle at first gsres@gsMarkerThicknessF = 1 gsres@gsMarkerSizeF = 0.01542 ; gsres@gsMarkerColor = "black" gsres@gsMarkerColor = "white" gsres@gsMarkerIndex = 16 ; Use filled dots for markers. ; if (j.eq.0) then dummy1 = gsn_add_polymarker(wks,contour_theta,x_lon,j_Altitude_GPS,gsres) txres = True txres@txFontHeightF = 0.01 txres@txBackgroundFillColor = "Transparent" interval = 1400/average_points text = gsn_add_text(wks,contour_theta,flt2string(hours_since_midnight_UTC_simple(::interval)),x_lon(::interval)+3, \ j_Altitude_GPS(::interval)+5 ,txres) dummy = new(num_distinct_markers,typeof(dummy1)) ; end if do j = 0, num_distinct_markers-1 if (.not.ismissing(y_new(j,0))) gsres@gsMarkerThicknessF = 1 gsres@gsMarkerColor = "black" gsres@gsMarkerIndex = 16 ; Use filled dots for markers. gsres@gsMarkerSizeF = 0.009087 gsres@gsMarkerColor = colors(j) ; if (j.eq.0) then ; dummy1 = gsn_add_polymarker(wks,contour_theta,x_new(j,:),y_new(j,:),gsres) ; printVarSummary(dummy1) ; end if dummy(j) = gsn_add_polymarker(wks,contour_theta,x_new(j,:),y_new(j,:),gsres) end if end do txres = True txres@txFontHeightF = 0.016881 txres@txBackgroundFillColor = "white" text_label = gsn_add_text(wks,contour_theta,"a",n_XMin+4, 10 ,txres) delete(y_new) delete(x_new) delete(colors) delete(levels) resvector = True resvector@vpHeightF = 0.253025 resvector@vpWidthF = 0.425 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 = 10. ; define vector ref mag resvector@vcRefLengthF = 0.025 ; define length of vec ref resvector@vcRefAnnoOrthogonalPosF = -1. ; move ref vector resvector@vcRefAnnoParallelPosF = 0.95 ; move ref vector resvector@vcMinDistanceF = 0.025 ; larger means sparser resvector@vcLineArrowHeadMaxSizeF = 0.0075 ; default: 0.05 (LineArrow), 0.012 (CurlyVector) resvector@gsnMaximize = False resvector@vcRefAnnoString2On = True; False resvector@vcRefAnnoString2 = "m s~S~-1" ; vector = wrf_vector(a,wks,va_plane,wa_plane*100,resvector) vector = wrf_vector(a,wks,ua_plane,wa_plane*100,resvector) vector2 = wrf_vector(a,wks,ua_plane,wa_plane*100,resvector) vector3 = wrf_vector(a,wks,ua_plane,wa_plane*100,resvector) vector4 = wrf_vector(a,wks,ua_plane,wa_plane*100,resvector) vector5 = wrf_vector(a,wks,ua_plane,wa_plane*100,resvector) vector6 = wrf_vector(a,wks,ua_plane,wa_plane*100,resvector) ;Contour terrain cross section ; print(opts_ter) opts_ter@trXMaxF = lon_plane(n_XMax) opts_ter@trXMinF = lon_plane(n_XMin) opts_ter@trXMaxF = n_XMax opts_ter@trXMinF = n_XMin opts_ter@trYMaxF = zcordinate(n_YMax)*1000 opts_ter@vpHeightF = 0.253025 opts_ter@vpWidthF = 0.425 opts_ter@xyLineThicknessF = 4.0 opts_ter@pmLegendSide = "Top" ; Change location of opts_ter@pmLegendParallelPosF = .01651084 ; move units right opts_ter@pmLegendOrthogonalPosF = -1.320386 opts_ter@pmLegendWidthF = 0.07081 ; Change width and opts_ter@pmLegendHeightF = 0.051 ; height of legend. opts_ter@pmLegendDisplayMode = "Always" ; create a legend opts_ter@lgPerimOn = False opts_ter@lgLabelFontHeightF = 0.0148 opts_ter@xyExplicitLegendLabels = (/"terrain","PBL top","lfc","lcl"/) opts_ter@xyLineColors = (/"blue","black","blue"/) opts_ter@xyDashPatterns = (/0, 12, 16/) opts_ter@gsnYRefLine = 0.0 ; opts_ter@gsnAboveYRefLineColor = "black" ter_plane_mean_plusref = new ((/2,dimsizes(ter_plane)/),float) ter_plane_mean_plusref(0,:) = (/ter_plane/) ter_plane_mean_plusref(1,:) = (/ter_plane+PBLH_plane/) ; contour_ter = gsn_csm_xy(wks,X_plane,ter_plane,opts_ter) ; contour_ter = gsn_csm_xy(wks,X_plane,ter_plane_mean_plusref,opts_ter) contour_ter = gsn_csm_y(wks,ter_plane_mean_plusref,opts_ter) opts_ter@pmLegendDisplayMode = "Never" ; create a legend contour_ter2 = gsn_csm_y(wks,ter_plane_mean_plusref,opts_ter) contour_ter3 = gsn_csm_y(wks,ter_plane_mean_plusref,opts_ter) contour_ter4 = gsn_csm_y(wks,ter_plane_mean_plusref,opts_ter) contour_ter5 = gsn_csm_y(wks,ter_plane_mean_plusref,opts_ter) contour_ter6 = gsn_csm_y(wks,ter_plane_mean_plusref,opts_ter) ; contour_ter = gsn_csm_xy(wks,X_plane,ter_plane,opts_ter) ; dummy1 = gsn_add_polymarker(wks,contour_ter,X_plane(id_d),0.08,gsres) ; j_Altitude_GPS = elevation_array(:,1)/(zmax*1000./dim(0)) ; x_lon = (lat-min(X_plane))/((max(X_plane)-min(X_plane))/dimsizes(X_plane)) ; print("max, min of X_plane"+ min(X_plane)+" "+max(X_plane)+"max, min of lat "+min(lat)+" "+max(lat) ) ; print("max, min of x_lon"+ min(x_lon)+" "+max(x_lon) ) ; dummy1 = gsn_add_polymarker(wks,contour_ter,lat,Altitude_GPS,gsres) ; MAKE PLOTS pltres = True pltres@NoTitles = True pltres@CommonTitles = True; False pltres@PanelPlot = True pltres@gsnPaperOrientation = "portrait" pltres@FramePlot = False pltres@PlotTitle = "Vertical cross-section along the C130 track" Times_char = a->Times CST_hr = mod(stringtoint(chartostring(Times_char(0,11:12)))-6+24, 24) if (ifiles.lt.6) then time_xhu = "Aug 3, "+ CST_hr+" CT" else time_xhu = "Aug 4, "+ CST_hr+" CT" end if ; plot = wrf_overlays(a,wks,(/contour_theta,vector/),pltres) plot = wrf_overlays(a,wks,(/contour_theta,vector,contour_ter/),pltres) setvalues contour_theta "vpHeightF" : 0.253025 "vpWidthF" : 0.425 "vpYF" : 0.9600 ; 0.9-ifile_loop*0.39 "vpXF" : 0.081 ;0.15;+ispecies*0.355 "tiMainString" : time_xhu ; "cross-section along the C130 track "+a->Times(0,6:12) + " UTC" "tiMainFontHeightF":0.0128 "tiMainOffsetYF" : -0.01 "gsnRightString" : "Vertical cross-section along the C130 track" ; res@TimeLabel ; "gsnPaperOrientation" : "portrait" end setvalues draw(contour_theta) opts_theta = opts_xy opts_theta@cnFillMode = "RasterFill" opts_theta@tmXBFormat = "f" opts_theta@pmLabelBarOrthogonalPosF = -0.07 opts_theta@cnFillOn = True opts_theta@gsnSpreadColors = True ; use full range of colormap opts_theta@lbTitleOn = True ; turn on title ; opts_theta@lbTitleString = " CO~B~2~N~~C~~V10~ ppmv" ; "O~B~3, ~N~ppbv" ; title string ; opts_theta@lbTitleString = " Temperature~C~~V10~~S~o~N~C" ; "O~B~3, ~N~ppbv" ; title string ; opts_theta@lbTitleString = " td~C~~V10~~S~o~N~C" ; "O~B~3, ~N~ppbv" ; title string ; opts_theta@lbTitleString = " theta~C~~V10~~S~o~N~C" ; "O~B~3, ~N~ppbv" ; title string opts_theta@lbTitlePosition = "Right" ; title position opts_theta@lbTitleFontHeightF= .0148 ; make title smaller opts_theta@lbLabelOffsetF = 0.01 opts_theta@lbTitleDirection = "Across" ; title direction ; opts_theta@lbTopMarginF = -0.051002 ; opts_theta@lbBottomMarginF = 0.21 opts_theta@lbTopMarginF = -0.151002 opts_theta@lbBottomMarginF = 0.321 opts_theta@lbLeftMarginF = 0.1251 opts_theta@lbRightMarginF = 0.1251 opts_theta@lbTitleOffsetF = 0.0004 opts_theta@lbTitleJust = "CenterCenter" opts_theta@tiYAxisString = " " ; Plotting options for Temperature opts_tc = opts_xy opts_tc@cnInfoLabelOrthogonalPosF = 0.00 opts_tc@ContourParameters = (/ 5. /) ; Get the contour info for the theta and temp opts_theta@gsnMaximize = False ; opts_theta@ContourParameters = (/ -24., 30., 2. /) ; tc_plane!0 = "height" ; tc_plane!1 = "lon" ; tc_plane&height = fspan(0,zmax,dim(0) ) ; tc_plane&lon = X_plane ; contour_tc = wrf_contour(a,wks,tc_plane,opts_theta) ; opts_theta@ContourParameters = (/ -34., 20., 2. /) ; td_plane!0 = "height" ; td_plane!1 = "lon" ; td_plane&height = fspan(0,zmax,dim(0) ) ; td_plane&lon = X_plane ; contour_tc = wrf_contour(a,wks,td_plane,opts_theta) thetae_plane!0 = "height" thetae_plane!1 = "lon" thetae_plane&height = fspan(0,zmax,dim(0) ) thetae_plane&lon = X_plane ; opts_theta@lbTitleString = "CO~B~2~N~Background~C~~V10~ppmv" ; "O~B~3, ~N~ppbv" ; title string ; contour_tc = wrf_contour(a,wks,CO2BCK_plane,opts_theta) theta_plane!0 = "height" theta_plane!1 = "lon" theta_plane&height = fspan(0,zmax,dim(0) ) theta_plane&lon = X_plane if ( WSP_crossSNpartial ) then ; = True opts_theta@ContourParameters = (/ 300., 325., 1. /) opts_theta@lbTitleString =" ~F5~q~F~~B~~N~, K"; " Thetae~C~~V10~~S~o~N~C" ; "O~B~3, ~N~ppbv" ; title string contour_tc = wrf_contour(a,wks,theta_plane,opts_theta) else opts_theta@lbTitleString =" ~F5~q~F~~B~e~N~, K"; " Thetae~C~~V10~~S~o~N~C" ; "O~B~3, ~N~ppbv" ; title string opts_theta@ContourParameters = (/ 320., 360., 2. /) contour_tc = wrf_contour(a,wks,thetae_plane,opts_theta) end if ; text_labelc = gsn_add_text(wks,contour_tc,"b",n_XMin+4, n_YMax-10 ,txres) text_labelc = gsn_add_text(wks,contour_tc,"b",n_XMin+4, 10 ,txres) gsres@gsMarkerThicknessF =6 gsres@gsMarkerSizeF = 0.007081002 gsres@gsMarkerColor = "blue" gsres@gsMarkerIndex = 7 ; Use filled dots for markers. ; dummy_front2 = gsn_add_polymarker(wks,contour_tc,237.5,6,gsres) txres = True txres@txFontHeightF = 0.012 txres@txBackgroundFillColor = "Transparent" ; text_front2 = gsn_add_text(wks,contour_tc,"front",237.5,.1 ,txres) ; plot2 = wrf_overlays(a,wks,(/contour_tc,vector2,contour_ter/),pltres) plot2 = wrf_overlays(a,wks,(/contour_tc,vector2,contour_ter2/),pltres) setvalues contour_tc "vpHeightF" : 0.253025 "vpWidthF" : 0.425 "vpYF" : 0.9600 ; 0.9-ifile_loop*0.39 "vpXF" : 0.55;+ispecies*0.355 "tiMainString" : time_xhu ;"cross-section along the C130 track "+a->Times(0,6:12) + " UTC" "tiMainFontHeightF":0.0128 "tiMainOffsetYF" : -0.01 end setvalues draw(contour_tc) ;;;;;;;;;;;;;;;;;add front symbol;;;;;;;;;;;;;;;;;; ; opts_theta@lbTitleString = "CO~B~2~N~Biogenic~C~~V10~ppmv" ; "O~B~3, ~N~ppbv" ; title string ; opts_theta@ContourParameters = (/ -6., 6., 1. /) ; opts_theta@tiYAxisString = "Height, km" ; contour_bio = wrf_contour(a,wks,CO2bio_plane,opts_theta) opts_theta@lbTitleString ="q, g kg~S~-1"; " QVAPOR~C~~V10~~S~o~N~C" ; "O~B~3, ~N~ppbv" ; title string opts_theta@ContourParameters = (/ 1, 15., .5 /) QVAPOR_plane!0 = "height" QVAPOR_plane!1 = "lon" QVAPOR_plane&height = fspan(0,zmax,dim(0) ) QVAPOR_plane&lon = X_plane contour_bio = wrf_contour(a,wks,QVAPOR_plane,opts_theta) txres@txFontHeightF = 0.016881 text_labelcc = gsn_add_text(wks,contour_bio,"c",n_XMin+4, 10 ,txres) plot3 = wrf_overlays(a,wks,(/contour_bio,vector3,contour_ter3/),pltres) setvalues contour_bio "vpHeightF" : 0.253025 "vpWidthF" : 0.425 "vpYF" : 0.640 ; 0.9-ifile_loop*0.39 "vpXF" : 0.081;+ispecies*0.355 "tiMainString" : "" "tiMainFontHeightF":0.0128 "tiMainOffsetYF" : -0.01 end setvalues draw(contour_bio) print("finsih bio") ;;;;;;;;;;;;;;;;;add front symbol;;;;;;;;;;;;;;;;;; ; opts_theta@lbTitleString = "CO~B~2~N~Anthropogenic~C~~V10~ppmv" ; "O~B~3, ~N~ppbv" ; title string ; opts_theta@ContourParameters = (/ -6., 6., 1. /) ; opts_theta@tiYAxisString = "" ; contour_Anth = wrf_contour(a,wks,CO2Anth_plane,opts_theta) opts_theta@lbTitleString ="WSP, m s~S~-1"; " WSP~C~~V10~~S~o~N~C" ; "O~B~3, ~N~ppbv" ; title string opts_theta@ContourParameters = (/ 2., 20., 1. /) WSP_plane!0 = "height" WSP_plane!1 = "lon" WSP_plane&height = fspan(0,zmax,dim(0) ) WSP_plane&lon = X_plane contour_Anth = wrf_contour(a,wks,WSP_plane,opts_theta) text_labeld = gsn_add_text(wks,contour_Anth,"d",n_XMin+4, 10 ,txres) plot4 = wrf_overlays(a,wks,(/contour_Anth,vector4,contour_ter4/),pltres) setvalues contour_Anth "vpHeightF" : 0.253025 "vpWidthF" : 0.425 "vpYF" : 0.640 ; 0.9-ifile_loop*0.39 "vpXF" : 0.55;+ispecies*0.355 "tiMainString" : "" "tiMainFontHeightF":0.0128 "tiMainOffsetYF" : -0.01 end setvalues draw(contour_Anth) print("finsih Anth") ;;;;;;;;;;;;;;;;;add front symbol;;;;;;;;;;;;;;;;;; ; opts_theta@lbTitleString = "CO~B~2~N~Resp~C~~V10~ppmv" ; "O~B~3, ~N~ppbv" ; title string ; opts_theta@ContourParameters = (/ -24., 24., 2. /) ; opts_theta@tiYAxisString = "Height, km" ; contour_Resp = wrf_contour(a,wks,CO2Resp_plane,opts_theta) gsn_define_colormap(wks,"BlAqGrWh2YeOrReVi22") ; select color map opts_theta@lbTitleString ="W, m s~S~-1"; " WSP~C~~V10~~S~o~N~C" ; "O~B~3, ~N~ppbv" ; title string opts_theta@ContourParameters = (/ -.5, .5, .05 /) wa_plane!0 = "height" wa_plane!1 = "lon" wa_plane&height = fspan(0,zmax,dim(0) ) wa_plane&lon = X_plane contour_Resp = wrf_contour(a,wks,wa_plane,opts_theta) text_labelf = gsn_add_text(wks,contour_Resp,"e",n_XMin+4, 10 ,txres) plot6 = wrf_overlays(a,wks,(/contour_Resp,vector5,contour_ter6/),pltres) setvalues contour_Resp "vpHeightF" : 0.253025 "vpWidthF" : 0.425 "vpYF" : 0.32 ; 0.9-ifile_loop*0.39 "vpXF" : 0.081 ;+ispecies*0.355 "tiMainString" : "" "tiMainFontHeightF":0.0128 "tiMainOffsetYF" : -0.01 end setvalues draw(contour_Resp) print("finsih Resp") ; opts_theta@lbTitleString = "CO~B~2~N~Photo~C~~V10~ppmv" ; "O~B~3, ~N~ppbv" ; title string ; opts_theta@ContourParameters = (/ -24., 24., 2. /) ; opts_theta@tiYAxisString = "" ; contour_Photo = wrf_contour(a,wks,CO2Photo_plane,opts_theta) gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ; select color map if (WSP_crossSNpartial) then opts_theta@lbTitleString ="QC, mg kg~S~-1"; " QVAPOR~C~~V10~~S~o~N~C" ; "O~B~3, ~N~ppbv" ; title string opts_theta@ContourParameters = (/ 2, 20, 2 /) QCLOUD_plane!0 = "height" QCLOUD_plane!1 = "lon" QCLOUD_plane&height = fspan(0,zmax,dim(0) ) QCLOUD_plane&lon = X_plane contour_Photo = wrf_contour(a,wks,QCLOUD_plane,opts_theta) else opts_theta@lbTitleString ="Kh, m~S~2~N~ s~S~-1"; " EXCH_H~C~~V10~~S~o~N~C" ; "O~B~3, ~N~ppbv" ; title string if(ifiles.lt.16) then opts_theta@ContourParameters = (/ .5, 10, .5 /) else opts_theta@ContourParameters = (/ 5, 100, 5 /) end if EXCH_H_plane!0 = "height" EXCH_H_plane!1 = "lon" EXCH_H_plane&height = fspan(0,zmax,dim(0) ) EXCH_H_plane&lon = X_plane contour_Photo = wrf_contour(a,wks,EXCH_H_plane,opts_theta) end if text_labele = gsn_add_text(wks,contour_Photo,"f",n_XMin+4, 10 ,txres) plot5 = wrf_overlays(a,wks,(/contour_Photo,vector6,contour_ter5/),pltres) setvalues contour_Photo "vpHeightF" : 0.253025 "vpWidthF" : 0.425 "vpYF" : 0.32 ; 0.9-ifile_loop*0.39 "vpXF" : 0.55 ; 081;+ispecies*0.355 "tiMainString" : "" "tiMainFontHeightF":0.0128 "tiMainOffsetYF" : -0.01 end setvalues draw(contour_Photo) print("finsih Photo") ; delete(lat_new) ; delete(lon_new) ; delete(colors) ; delete(levels) delete(res) end do ; END OF TIME LOOP delete(ter_plane) delete(wa_plane) delete(ua_plane) delete(WSP_plane) delete(QVAPOR_plane) delete(EXCH_H_plane) delete(tc_plane) delete(theta_plane) delete(va_plane) delete(CO2total_plane); = wrf_user_intrp3d(CO2total,z,"v",plane,angle,opts) delete(zcordinate) frame(wks) ; system("convert -density 300 -rotate -90 -resize 1000x1000 -trim "+figurename+".eps /nsftor/xhu/public_html/HurricaneImpactonO3/WRFV3.6.1/YSU/wrfNARR3dWSM6_CONUS_UCM_YSU.2011080700/"+figurename+".png") ; system("convert -density 300 -resize 1000x1000 -trim "+figurename+".eps /nsftor/xhu/public_html/HurricaneImpactonO3/WRFV3.6.1/YSU/wrfNARR3dWSM6_CONUS_UCM_YSU.2011080700/"+figurename+".png") print("finish "+figurename+".eps") print("finish working on " + files(ifiles)) end do ; ifiles end