;Written by Xiaoming Hu,CAPS, OU, Oct. 3, 2011 ; 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" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" begin idomain = 1 S54511_Lat = 39.8 S54511_Lon = 116.4667 latitude_3S54511 = (/S54511_Lat, S54511_Lat, S54511_Lat/) longitude_3S54511 = (/S54511_Lon, S54511_Lon, S54511_Lon/) names_3S54511 = (/"S54511", "S54511", "S54511"/) lat = S54511_Lat lon = S54511_Lon timewithsounding = (/5, 11, 23, 29/) ; soundingfiles = systemfunc("ls Z_UPAR_I_59981_201?????05*.SECOND") ; soundingfiles = systemfunc("ls Z_UPAR_I_59981_201?????23*.SECOND") soundingfiles = systemfunc("ls Z_UPAR_I_59981_201?????11*.SECOND") ; also turn off adding the PBL line wrfoutfiles = systemfunc("ls /oasis/scratch/comet/xhu2/temp_project/Run/CO2_and_otherGHG/WRFV3.9.1.1/China/wrfchem3.9.1.1_R2_China2Shanghai_10mb_Hil3ReShrubRES_addOce_restoreDF_ODIAC_CT2017.2016010100/wrfout_d0"+idomain+"*2016-12*12:00:00") resp = True ; polyline mods desired resp@gsLineColor = "red" ; color of lines resp@gsLineThicknessF = 4.0 ; thickness of lines resp@gsLineDashPattern= 2 ; thickness of lines tres = True tres@txFontHeightF = 0.01322 fmap = addfile(wrfoutfiles(0)+".nc","r") lat2d=fmap->XLAT(0,:,:) ;(Time, south_north, west_east) lon2d=fmap->XLONG(0,:,:) HGT= fmap->HGT(0,:,:);TerrainHeight ZNU=fmap->ZNU(0,:) nlayers = dimsizes(ZNU) data_latloncontain = asciiread(soundingfiles(0),-1,"float"); (soundingfiles(ihours), 12, "float", 7) lat = data_latloncontain(3) lon = data_latloncontain(2) JI_U = getind_latlon2d(lat2d,lon2d, lat,lon) dimdomain = dimsizes(HGT) upbound = dimdomain(0)-2 rightbound = dimdomain(1)-2 print("rightbound"+rightbound+" upbound"+upbound) if (.not.all(ismissing(soundingfiles)).and.JI_U(0,0).gt.2.and.JI_U(0,0).lt.upbound.and.JI_U(0,1).gt.0.and.JI_U(0,1).lt.rightbound) then print("find lat, lon"+lat+" "+lon) do ihours = 0, dimsizes(soundingfiles)-2 ; figurename = "wrfout_d0"+idomain+"_profile_Theta_q_5z_Yongxingdao_"+ihours ; figurename = "wrfout_d0"+idomain+"_profile_Theta_q_23z_"+ihours figurename = "wrfout_d0"+idomain+"_profile_Theta_q_11z_Yongxingdao_"+ihours ; figurename = "wrfout_d0"+idomain+"_profile_Theta_withPBL_q_"+ihours f=addfile(wrfoutfiles(ihours)+".nc","r") Times_char = f->Times Times_char(0,13) = Times_char(0,4) Times_char(0,16) = Times_char(0,4) ua = wrf_user_getvar(f,"ua",0) ; u on mass points va = wrf_user_getvar(f,"va",0) ; v on mass points WSP = sqrt(ua^2+va^2) height = wrf_user_getvar(f,"height",0) Jmin_u = JI_U(0,0)-1 Jmax_u = JI_U(0,0)+1 Imin_u = JI_U(0,1)-1 Imax_u = JI_U(0,1)+1 ; print("Terrain height="+dim_avg_Wrap(dim_avg_Wrap(HGT(Jmin_u:Jmax_u,Imin_u:Imax_u))) ) Height_UrbanAve = dim_avg_Wrap(dim_avg_Wrap(height(:,Jmin_u:Jmax_u,Imin_u:Imax_u))) - dim_avg_Wrap(dim_avg_Wrap(HGT(Jmin_u:Jmax_u,Imin_u:Imax_u))) WSP_urbanArea = WSP(:, JI_U(0,0), JI_U(0,1)) Tpotential_urbanArea = f->T(0,:, JI_U(0,0), JI_U(0,1)) QVAPOR_urbanArea = f->QVAPOR(0,:, JI_U(0,0), JI_U(0,1)) a4_height = 29.7 ; in centimeters, if my a4_width = 23.0 ; reference is correct cm_per_inch = 2.54 res1 = True ; res1@gsnFrame = False ; don't draw yet ; res1@mpShapeMode = "FreeAspect" res1@vpHeightF =0.42 res1@vpWidthF = 0.42 ; res1@vpYF =0.65 ; res1@gsnPaperHeight = a4_height/cm_per_inch ; res1@gsnPaperWidth = a4_width/cm_per_inch res1@gsnPaperOrientation = "portrait" ; res1@gsnMaximize = True res1@gsnPaperMargin = 0.01 ; res1@tiXAxisString = "~F8~q~F21~~B~v~N~, K";"O~B~3~N~, ppb";"Potential Temperature, K"; ; res1@xyExplicitLegendLabels = leg ; explicit labels created above res1@xyDashPatterns = (/0,0,16,14/) ; 0 is default (solid) res1@xyLineColors = (/"red","red","black","red","Magenta2","red","brown"/) res1@xyMarkLineModes = (/"Markers","Lines","Lines","Lines","Lines" /) res1@xyMarkers = (/16,16/) ; choose type of marker ; res1@xyMarkerColor = "Turquoise3" res1@xyMarkerColor = "Cyan3" res1@xyMarkerColor = "red" res1@xyMarkerSizeF = 0.00855 res1@xyLineThicknessF = 8 res1@pmLegendSide = "Top" ; Change location of res1@pmLegendParallelPosF = .21983 ; move units right res1@pmLegendOrthogonalPosF = -.36291 res1@pmLegendWidthF = 0.08102 ; Change width and res1@pmLegendHeightF = 0.135 ; height of legend. res1@tiXAxisFuncCode = "~" ; res1@tmXLMode = "Manual" ; res1@trXMaxF = Tpotential_ea_sitesAve(1, 20) res1@trYMaxF = 3. if (ihours.eq.29) then res1@trXMaxF = 310 ; res1@trXMaxF = 330 res1@trXMinF = 301 end if res1@tmYLValues = fspan(0,5,11) res1@tmYLMinorValues= fspan(0,5,21) res1@tmYRMinorValues= fspan(0,5,21) res1@tmYLLabels = fspan(0,5,11) res1@tmYLMode = "Explicit" res1@trYMinF = 0 res1@tmYLLabelStride = 1 res1@tmYRLabelStride = 2 ; res1@tmXBLabelStride = 2 res1@tmXTLabelStride = 2 res1@tmYLFormat = "f" res1@tmYRFormat = "f" res1@tmXBFormat = "f" res1@tmYLTickSpacingF = 0.5 res1@tmXBTickSpacingF = 5.0 res1@tmYRMajorOutwardLengthF = 0 res1@tmYLMajorOutwardLengthF = 0 res1@tmXBMajorOutwardLengthF = 0 res1@tmXBMinorOutwardLengthF = 0 res1@tmYLMinorOutwardLengthF = 0 res1@tmXTMinorOutwardLengthF = 0 res1@tmYRMinorOutwardLengthF = 0 res1@tmYLMinorOn = True res1@tmYRMinorOn = True res1@tmXBMinorOn = True res1@tmXTMinorOn = True res1@lgPerimOn=False res1@lgLabelFontHeightF = 0.0205 res1@tiXAxisFontHeightF= 0.028 res1@tiYAxisFontHeightF= 0.028 res1@tmXBLabelFontHeightF = 0.028 res1@tmYLLabelFontHeightF = 0.028 res1@lgLabelFontHeightF = 0.0205 res1@tiXAxisFontHeightF= 0.0158 res1@tiMainFontHeightF = 0.0158 res1@tiYAxisFontHeightF= 0.0158 res1@tmXBLabelFontHeightF = 0.0158 res1@tmYLLabelFontHeightF = 0.0158 res1@tmXMajorGrid = True ; implement x grid res1@tmXMajorGridThicknessF = 1.0 ; 2.0 is default res1@tmXMajorGridLineDashPattern = 2 ; select short dash lines res1@pmTickMarkDisplayMode = "Always" ; turn on tickmarks ; res1@tmXTOn = False ; turn off top labels ; res1@tmYROn = False ; turn off right labels ; res1@trYReverse = True res1@pmLegendDisplayMode = "Never" ; create a legend ; res1@xyExplicitLegendLabels = (/"05 UTC","09 UTC"/) ; res1@lgItemOrder = (/3,2,1,0/) res2=res1 res2@xyMarkLineModes = (/"MarkLines","Lines","Lines","Lines","Lines" /) res2@xyMarkerColor = "black" res2@xyLineColors = (/"black","red","black","red","Magenta2","red","brown"/) res2@xyDashPatterns = (/16,0,16,14/) ; 0 is default (solid) ; res2@tmXTMode = "Explicit" ; res2@tmXTValues = ispan(0,20,2) ; res2@tmXTLabels = ispan(0,20,2) res2@tmYRLabelStride = 2 res2@tmYRTickSpacingF = 0.5 res2@tiXAxisString = "" res2@tiYAxisString = "" res2@tmYRMode = "Manual" wks_type = "png" wks_type@wkWidth = 811 wks_type@wkHeight = 811 ; wks_type@wkWidth = 1011 ; wks_type@wkHeight = 1011 wks = gsn_open_wks(wks_type,figurename) ; wks = gsn_open_wks("eps" ,figurename) ; eps,pdf,x11,ncgm,eeps res2@tiDeltaF = 0.5 ; res1@tiMainOffsetYF=-0.01 res2@tmYRLabelsOn = False res2@tmXTFormat = "f" ; print(Tpotential_ea_sitesAve(:,:3)) ; print("gonna reading "+soundingfiles(ihours)) res2@xyExplicitLegendLabels = (/"WRF D"+idomain/) data_second = readAsciiTable(soundingfiles(ihours), 12, "float", 7) file_char = stringtocharacter(soundingfiles(ihours)) ; print(file_char) res1@xyExplicitLegendLabels = (/"Radiosonde"/); chartostring(file_char(9:13)) data_second@_FillValue = -999. ; Height_second = sin(data_second(:,4)*3.1415926535/180.)*data_second(:,6) data_all = asciiread(soundingfiles(ihours), -1, "float") height_site = data_all(4) delete(data_all) Height_second = (data_second(:,11)-height_site)/1000. index_bad = ind(Height_second.eq.0) ; print(Height_second(index_bad)) if(.not.all(ismissing(index_bad))) Height_second(index_bad) = Height_second@_FillValue end if delete(index_bad) Press = data_second(:,2) Press = where(Press.ne.0, Press, Press@_FillValue) Theta_second = (data_second(:,1)+273.15)*(1000./Press)^0.286 ; print(Theta_second(:10)) ; print(data_second(:10,1)) ;;;;;;;;;;;;;;;;;;;;diagnose PBLH from LBand Radiosound;;;;;;;;;;;;;;;;;;;;;;;;;; minPT_Obs = min(Theta_second) ilevel = minind(Theta_second)+1 ismissing_xhu = all(ismissing(Theta_second(ilevel))) do while (ismissing_xhu) ; print("ilevel "+ilevel +" "+Theta_second(ilevel)) ilevel = ilevel + 1 ismissing_xhu = all(ismissing(Theta_second(ilevel))) end do ; ; if (icases.eq.15) then ; print("ilevel "+ilevel +" "+Theta_second(ilevel)) ; print("(Theta_second(ilevel)-minPT_Obs)=" + (Theta_second(ilevel)-minPT_Obs)) ; end if do while ((Theta_second(ilevel)-minPT_Obs).lt.1.5) ilevel = ilevel + 1 ismissing_xhu = all(ismissing(Theta_second(ilevel))) do while (ismissing_xhu) ; print("ilevel "+ilevel +" "+Theta_second(ilevel)) ilevel = ilevel + 1 ismissing_xhu = all(ismissing(Theta_second(ilevel))) end do ; end do PBLH_1_5_Obs = Height_second(ilevel-1) ;+ (1.5-(Theta_second(ilevel-1)-minPT_Obs))/(Theta_second(ilevel)-Theta_second(ilevel-1))* (Height_second(ilevel)-Height_second(ilevel-1)) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; q_sounding = mixhum_ptrh (Press, data_second(:,1)+273.15, data_second(:,3),-1) ; print(data_second(:10,3)) res1@trXMinF = 290.; max(Tpotential_ea_sitesAve(1, :29)) res1@trXMaxF = 330. ; min((/Tpotential_ea_sitesAve(1, 0)-6,TH2_Mesonet_mean-0.5/)) res1@trYMaxF = 5. res2@trYMaxF = 5. res1@trXMinF = 258.; for Shenyang winter max(Tpotential_ea_sitesAve(1, :29)) res1@trXMaxF = 312. ; min((/Tpotential_ea_sitesAve(1, 0)-6,TH2_Mesonet_mean-0.5/)) res2@trXMinF = 258.; for Shenyang winter max(Tpotential_ea_sitesAve(1, :29)) res2@trXMaxF = 312. ; min((/Tpotential_ea_sitesAve(1, 0)-6,TH2_Mesonet_mean-0.5/)) ; res2@trXMaxF = max(Tpotential_ea_sitesAve(1, :29)) ; res2@trXMinF = min((/Tpotential_ea_sitesAve(1, 0)-6,TH2_Mesonet_mean-0.5/)) res1@vpXF = 0.1 res1@tiXAxisString = "Potential temperature, K";"Potential Temperature, K"; res1@tiYAxisString = "Height, km"; "PM~B~2.5~N~, ~F33~m~F~g m~S~-3" ; label bottom axis with units res1@pmLegendDisplayMode = "Always" ; create a legend res2@pmLegendDisplayMode = "Always" ; create a legend res1@tiMainString = chartostring(file_char(15:26)) +" at " + chartostring(file_char(9:13)) + " Yongxingdao_Hainan" if (file_char(9:13).eq."55299") then res1@trXMinF = 310.; max(Tpotential_ea_sitesAve(1, :29)) res1@trXMaxF = 345. ; min((/Tpotential_ea_sitesAve(1, 0)-6,TH2_Mesonet_mean-0.5/)) end if ; plot = gsn_csm_xy(wks,Theta_second,Height_second,res1) plot = gsn_csm_x2y2(wks,Theta_second,Tpotential_urbanArea+300,Height_second,Height_UrbanAve/1000, res1, res2) resp@gsLineColor = "red" ; color of lines ; deltaPBLH = min((/abs(PBLH-PBLH_1_5_Obs), abs(PBLH_1_5-PBLH_1_5_Obs)/)) res1@vpXF = 0.53 res1@trXMaxF = 20. res2@trXMaxF = 20. res1@trXMaxF = 8.; for Shenyang winter res2@trXMaxF = 8. res1@trXMaxF = 12.; for Shenyang winter Oct res2@trXMaxF = 12. res1@trXMinF = .02 res2@trXMinF = .02 res1@tmYLLabelsOn = False res2@tmXTLabelsOn = True res1@tiXAxisString = "Mixing ratio, g kg~S~-1" res1@tiYAxisString = ""; "PM~B~2.5~N~, ~F33~m~F~g m~S~-3" ; label bottom axis with units res1@pmLegendDisplayMode = "Never" ; create a legend res2@pmLegendDisplayMode = "Never" ; create a legend ; print(q_sounding(:)) ; print(Height_second(:)) ; print(Height_second) ; print(q_sounding) ; res1@xyMarkerColor = "blue" res1@tiMainString =chartostring(f->Times(0, :12)) + "UTC d0"+idomain ; chartostring(file_char(15:24)) +" at " + chartostring(file_char(9:13)) print("working on "+ chartostring(f->Times(0, :12))+ " domain"+idomain +" Yongxingdao_Hainan") plot = gsn_csm_x2y2(wks,q_sounding,QVAPOR_urbanArea*1000, Height_second,Height_UrbanAve/1000, res1, res2) res1@tmYLLabelsOn = True delete(data_second) delete(Height_second) delete(Press) delete(Theta_second) delete(q_sounding) gsres = True gsres@gsMarkerThicknessF = 1 gsres@gsMarkerIndex = 13 gsres@gsMarkerSizeF = 0.015 gsres@gsMarkerColor = "pink" ; print(TH2_Mesonet_mean) frame(wks) ; system("convert -trim -resize 100% "+figurename+".eps /nsftor/xhu/public_html/BeijingO3_MaZhiQiang/Simulations/3.6/wrfYMiao_YSU.20100912/"+figurename+".png") ; system("rm "+figurename+".eps") ; print("finish "+figurename+".eps") end do ; ihour end if ; if (.not.all(ismissing(soundingfiles)).and.JI_U(0,0).gt.2.and.JI_U(0,0).lt.264.and.JI_U(0,1).gt.0.and.JI_U(0,1).lt.440) then end