;Written by Xiaoming Hu, Penn State, Jan, 2010 ; contact yuanfangcan@gmail.com load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl" function theta_eqv(T:numeric, X:numeric, p0:numeric, p:numeric) ; ; calculate equivalent potential temperature ; ; T : Temperature [K] (time,lev,lat,lon) ; X : mixing ratio of water vapor in air [kg/kg] (time,lev,lat,lon) ; p0 : reference pressure [ usually 1000 mb or 100000 Pa ] ; p ; pressure at each level [must be same units as p0 ] local T_e, L_v, c_pd, R, Rcpd, dimT, dimX, P begin L_v = 2400. ; latent heat of evaporation [approx 2400 kJ/kg at 25C] c_pd = 1004. ; specific heat at constant pressure for air [approx 1004 J/(kg-K)] R = 287. ; specific gas constant for air [J/(kg-K)] Rcpd = R/c_pd dimT = dimsizes(T) dimX = dimsizes(X) if (dimsizes(dimT) .ne. dimsizes(dimX)) then print("theta_eqv: rank of T .ne. rank of X") exit end if ; P = conform(T,p,1) ; make p same shape/rank/size as T P = p T_e = T + (L_v/c_pd)*X ; common approximation theta_e = T_e*(p0/P)^Rcpd theta_e_at_long_name = "equivalent potential temperature" theta_e_at_units = "K" copy_VarCoords(T, theta_e) ; assign coordinates return(theta_e) end begin sites = (/"TMFA","TMFB","WHFA","GRFA","BHFA"/) Colors = (/"red", "green","blue","Turquoise3", "Cyan3","Magenta2","Magenta2","green","green"/) variables=(/"IL","ustar","Qh","QhTc","Qe","FCo2","U","V","W","T","Tc","H2O","CO2","WSPD","WDIR","UU","VV","WW","TT","TCTC","HH","CC","wT","wTC","wq","wC","rho_a","rho_v","cp","ea","es","T2","RH2","P","Vane","Kdn","Kup","Ldn","Lup","Tkz","Qstar","QG","Tsoil","Moist","SW","Rain"/) variables_show = variables variables_unit=(/"n","m s~S~-1~N~","W/m^2","W/m^2","W/m^2","umol/m^2/s","m/s","m/s","m/s","C","C","g/m^3","mmol/m^3","m/s","Deg","^2/100","^2/100","^2/100","^2/100","^2/100","^2/100","^2/100","cm*T/s","cm*T/s","cm*q/s","cm*C/s","kg/m^3","g/m^3","J/K/kg","hPa","hPa","C","%","hPa","deg","W/m^2","W/m^2","W/m^2","W/m^2","C","W/m^2","W/m^2","C","0/0","N/A","mm"/) do isites=4,4 ; 3, dimsizes(sites)-1 data_array = readAsciiTable("indianau_energy"+sites(isites)+"-0317800-0321324.dat",50,"double",76) data_array@_FillValue = -999. print(data_array(0,:)) print(data_array(:10,0)) ; do ivar = 0, dimsizes(variables)-1 ; ; do ivar = 1,1 do ivar = 13,13; 11,11 do iday = 181,181; 212; 239 iday_number = iday-181 ; figurename = "TS_"+variables(ivar)+"_at"+sites(isites)+"_"+iday_number figurename = "TS_WindVector2days_at"+sites(isites)+"_"+iday_number Day_Hour_min_sec = data_array(:,0)+data_array(:,1)/24.0+data_array(:,2)/24.0/60+data_array(:,3)/24.0/60./60. - 2003000. -6./24. ; convert to CST ; print(Day_Hour_min_sec) Data_to_plot = data_array(:,ivar+4) WSP = data_array(:,ivar+4) WDR = data_array(:,ivar+5) pi = 3.1415926535 u =WSP*(- sin (WDR*pi/180)) v =WSP*(- cos(WDR*pi/180)) ; index_weird = ind(Data_to_plot.lt.5) ; print(index_weird) ; print(data_array(694,:)) dims = dimsizes(data_array) a4_height = 29.7 ; in centimeters, if my a4_width = 23.0 ; reference is correct cm_per_inch = 2.54 res1 = True ; plot mods desired if (variables(ivar).eq."ustar") then res1@trYMinF =0 res1@trYMaxF =1.2 end if if (variables(ivar).eq."NO") then res1@trYMinF =0 res1@trYMaxF =25 end if if (variables(ivar).eq."NO2") then res1@trYMinF = 0 res1@trYMaxF =35 end if if (variables(ivar).eq."NOx") then res1@trYMinF = 0 res1@trYMaxF =55 end if if (variables(ivar).eq."Ozone") then res1@trYMinF =0 res1@trYMaxF =120 end if if (variables(ivar).eq."WDIR") then res1@trYMinF =0 res1@trYMaxF =361 res1@tmYLMode = "Explicit" res1@tmYLValues = ispan(0,360,90) res1@tmYLLabels = ispan(0,360,90) res1@tmYMajorGrid = True ; implement x grid res1@tmYMajorGridThicknessF = 0.8 ; 2.0 is default res1@tmYMajorGridLineDashPattern = 2 ; select short dash lines end if if (variables(ivar).eq."SO2") then res1@trYMinF =-0.6 res1@trYMaxF = 12 end if if (variables(ivar).eq."CO") then res1@trYMinF =0 if (iday.eq.236) then res1@trYMaxF =1.8 else res1@trYMaxF =0.5 end if end if res1@gsnFrame = False ; don't draw yet ; res1@gsnDraw = False ; don't advance frame ; res1@mpShapeMode = "FreeAspect" res1@vpHeightF =0.3 res1@vpWidthF = 0.8 res1@vpXF = 0.15 res1@tiYAxisOffsetXF = 0.02 ; res1@vpYF =-0.9 ; res1@gsnPaperOrientation = "landscape" ; res1@gsnMaximize = True if (iday.eq.181) then res1@trXMinF = 182. res1@trXMaxF = 213. index_to_plot = ind(Day_Hour_min_sec.le.213.and.Day_Hour_min_sec.ge.182) else res1@trXMinF = iday ;223 res1@trXMaxF = iday+2;240 index_to_plot = ind(Day_Hour_min_sec.le.(iday+2).and.Day_Hour_min_sec.gt.iday) ; res1@trXMinF = iday+0.5 ;223 ; res1@trXMaxF = iday+2.5;240 end if res1@gsnPaperMargin = 0.01 res1@tiYAxisString = "Wind vectors"; variables_show(ivar) + ", "+variables_unit(ivar) res1@tmYLLabelsOn = False res1@tmYLMajorOn = False res1@xyDashPatterns = (/0,16,0,16,0,4/) ; 0 is default (solid) res1@xyMarkers = (/16,16,16,16,16,16/) ; choose type of marker ; res1@xyMarkerColor = "Turquoise3" res1@xyLineColors = (/"red","black","green","blue","Purple4","Turquoise3","Purple4","Magenta2","red","brown"/) res1@xyMarkerColors = Colors(0) res1@xyMarkerSizeF = 0.018 res1@xyLineThicknessF = 5 res1@pmLegendSide = "Top" ; Change location of res1@pmLegendParallelPosF = .68 ; move units right res1@pmLegendWidthF = 0.4 ; Change width and res1@pmLegendHeightF = 0.18 ; height of legend. res1@tiXAxisFuncCode = "~" ; res1@tmYLMode = "Manual" ; res1@trYMinF = 294.0 ; res1@tmYLTickStartF= 294.0 ; res1@tmYLTickEndF = 304.0 res1@tmYLTickSpacingF = 2.0 res1@pmTickMarkDisplayMode = "Always" ; turn on tickmarks res1@tmXTOn = False ; turn off top labels res1@tmYROn = False ; turn off right labels res1@lgPerimOn = False res1@lgLabelFontHeightF = 0.0205 res1@tiXAxisFontHeightF= 0.028 res1@tiYAxisFontHeightF= 0.028 res1@tmXBLabelFontHeightF = 0.028 res1@tmYLLabelFontHeightF = 0.028 res1@tmYLFormat = "f" res1@pmLegendSide = "Top" ; Change location of res1@pmLegendParallelPosF = .85 ; move units right res1@pmLegendOrthogonalPosF = -1.15 ; move units down res1@pmLegendWidthF = 0.1 ; Change width and res1@pmLegendHeightF = 0.14 ; height of legend. res1@tmXBMode = "Explicit" if (iday.eq.181) then res1@tmXBValues = ispan(182,212,1) res1@tmXBLabels = "7/"+ispan (1,31,1) res1@tiXAxisString = "Date";;"DOY"; "PM~B~2.5~N~, ~F33~m~F~g m~S~-3" ; label bottom axis with units res1@tmXBLabelStride = 4 res1@xyMarkLineModes = (/"Lines","MarkLines","MarkLines","MarkLines","MarkLines","MarkLines" /) else res1@xyMarkLineModes = (/"Lines","Markers","Markers","Markers","Markers","Markers" /) res1@tiXAxisString = "Time, CST";;"DOY"; "PM~B~2.5~N~, ~F33~m~F~g m~S~-3" ; label bottom axis with units res1@tmXBLabelStride = 4 res1@tmXBValues = fspan(182,213,1+(213-182)*8) res1@tmXBLabels = (/"7/1","","6","","12","","18","","7/2","","6","","12","","18","","7/3","","6","","12","","18","","7/4","","6","","12","","18","","7/5","","6","","12","","18","","7/6","","6","","12","","18","","7/7","","6","","12","","18","","7/8","","6","","12","","18","","7/9","","6","","12","","18","","7/10","","6","","12","","18","","7/11","","6","","12","","18","","7/12","","6","","12","","18","","7/13","","6","","12","","18","","7/14","","6","","12","","18","","7/15","","6","","12","","18","","7/16","","6","","12","","18","","7/17","","6","","12","","18","","7/18","","6","","12","","18","","7/19","","6","","12","","18","","7/20","","6","","12","","18","","7/21","","6","","12","","18","","7/22","","6","","12","","18","","7/23","","6","","12","","18","","7/24","","6","","12","","18","","7/25","","6","","12","","18","","7/26","","6","","12","","18","","7/27","","6","","12","","18","","7/28","","6","","12","","18","","7/29","","6","","12","","18","","7/30","","6","","12","","18","","7/31","","6","","12","","18","","8/1"/) end if res1@tmYLLabelStride = 2 res1@tmYLOn = False res1@tmYLMinorOn = False ; res1@pmLegendDisplayMode = "Always" ; create a legend res1@xyExplicitLegendLabels = (/""/) wks = gsn_open_wks("eps" ,figurename) ; eps,pdf,x11,ncgm,eeps gsn_define_colormap(wks,"BlAqGrYeOrRevi200") res1@gsnLeftString = "" ;files(ifiles) res1@gsnRightString = "IU site: "+sites(isites) ; "" ; "1999 Kwajalein" plot = gsn_csm_xy (wks,Day_Hour_min_sec,Day_Hour_min_sec*0,res1) ; create plot wmsetp("vrs - reference vector size",12.9) wmsetp("vrn - NDC size of a reference vector size",0.3) wmsetp("vcc - vector color",1) wmsetp("vcw - vector linewidth scale factor",4.) y = Day_Hour_min_sec y = (/0/) wmvect(wks, dble2flt(Day_Hour_min_sec(index_to_plot)), dble2flt(y(index_to_plot)), dble2flt(u(index_to_plot)), dble2flt(v(index_to_plot)) ) wmvlbl(wks, 0.93 ,0.51) ; plot = gsn_csm_xy (wks,Day_Hour_min_sec,Data_to_plot,res1) ; create plot ; getvalues plot ; retrieve some of the plot resources ; "tmXBValues" : tmXBValues ; values used by NCL at major tick marks ; end getvalues ; print(tmXBValues) ; tmXBLabels = sprintf("%5.1f",res1@tmXBLabels ) ; setvalues plot ; "tmXBMode" : "Explicit" ; "tmXBValues" : tmXBValues ; values used by NCL at major tick marks ; "tmXBLabels" : tmXBLabels ; end setvalues ;variables = (/"MeanWSP", "FrictionVelocity", "sigma_w", "turbulence_intensity", "VirtualheatFlux", "MomentumFlux"/) ; draw(plot) frame(wks) system("convert -trim -resize 100% "+figurename+".eps /nsftor/xhu/public_html/WRF-UCM/EnergyBudget/"+figurename+".png") system("rm "+figurename+".eps") print("finish plotting "+figurename) ; delete(tmXBValues) ; delete(tmXBLabels) delete(res1@trYMinF) delete(res1@trYMaxF) delete(res1@tmXBValues) delete(res1@tmXBLabels) delete(index_to_plot) end do ; iday end do ; ivar delete(data_array) end do ; isites end