; 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 year = 2011 do imonth = 8,8; 1,12; 6; 7,12; 1,10 ; 8,9 ; 10,10 ; 1, 12 if(imonth.lt.10) then folder = year+"0"+imonth ; ARMfiles = systemfunc("ls sgpdwlwind.c1."+year+"0"+imonth+"*.cdf") ARMfiles = systemfunc("ls sgpLidarwind.c1."+year+""+"0"+imonth+"*.cdf") else folder = year+""+imonth ; ARMfiles = systemfunc("ls sgpdwlwind.c1."+year+""+imonth+"*.cdf") ARMfiles = systemfunc("ls "+folder+"/sgpLidarwind.c1."+year+""+imonth+"*.cdf") end if system("mkdir -p /nsftor/xhu/public_html/ARM/915rwpwindconC1/"+folder) species_infile = (/ "WSPD","WDIR"/) do ipower = 0,0 ; 3 ;The lower power modes are more used for the boundary layer, while higher power are used for further up in the atmosphere. To basically change the power, the pulse length is increased (thus resolution decreases and the height of the first range gate increases). do ifiles = 0,0; 0, dimsizes(ARMfiles)-2 file_to_open = ARMfiles(ifiles) file_to_open2 = ARMfiles(ifiles+1) print("will work on "+file_to_open) print("will work on "+file_to_open2) f=addfile(file_to_open+".nc","r") f2=addfile(file_to_open2+".nc","r") species = (/ "w_spd","w_dir"/) ; Real_DOY = mod(DOY, 1000) do ispecies = 0 , 0 ; dimsizes(species)-1 figurename = "ARMlidartime_height_"+species_infile(ispecies)+"_nightCenter_"+ifiles wks = gsn_open_wks("eps" ,figurename) ; gsn_define_colormap(wks, "BlAqGrYeOrRevi200") time_p1=f->time time_p2= cd_convert( f2->time, time_p1@units ) ; print(time_p1) ; print((f2->time)/3600.) ; print(time_p2) time = array_append_record (time_p1, time_p2, 0) delete(time_p1) delete(time_p2) height=f->height species_value_p1 = f->$species(ispecies)$ species_value_p2 = f2->$species(ispecies)$ species_value = transpose(array_append_record (transpose(species_value_p1), transpose(species_value_p2), 0)) delete(species_value_p1) delete(species_value_p2) species_value&height = height/1000. ; the 2nd dimension time_inhour=time/3600. ; already in CST; - 6. ; convert to CST print(time_inhour) species_value&time = time_inhour ; the 1st dimen delete(time_inhour) dims=dimsizes(species_value) species_value_1d = ndtooned(species_value) ; index_missing = ind(abs(species_value_1d).gt.100) index_missing = ind(abs(species_value_1d).eq."nan") if (.not.all(ismissing(index_missing))) then ; print(species_value_1d(index_missing)) species_value_1d(index_missing) = species_value_1d@_FillValue ; print(species_value_1d(index_missing)) end if delete(index_missing) species_value = (/onedtond(species_value_1d,(/dims(0),dims(1)/))/) delete(species_value_1d) res1 = True ; plot mods desired res1@gsnSpreadColors = True res1@vpHeightF =0.2 ; res1@vpWidthF = 0.92 res1@vpWidthF = 0.6 res1@vpXF = 0.2 res1@tiYAxisOffsetXF = 0.015 res1@tmYLLabelDeltaF = -0.99 res1@gsnMaximize = True res1@gsnPaperOrientation = "portrait" res1@trYMaxF = 3. res1@gsnPaperMargin = 0.01 res1@tiYAxisString = "Height, km";+temp@units; "PM~B~2.5~N~, ~F33~m~F~g m~S~-3" ; label bottom axis with units res1@tiXAxisString = "Time (UTC)";" ;print(res1@sfXArray) res1@pmLegendSide = "Top" ; Change location of res1@pmLegendParallelPosF = .8 ; move units right res1@pmLegendOrthogonalPosF = -0.5 ; move units down res1@pmLegendWidthF = 0.25 ; Change width and res1@pmLegendHeightF = 0.1 ; height of legend. res1@tiXAxisFuncCode = "~" ; if(iDOY.eq.177) then ; res1@trXMinF = 182. ; res1@trXMaxF = 213. unit=time@units unit_char = stringtochar(unit) res1@tiXAxisString = "CST";"CST Hour "+ chartostring(unit_char(8:)) ; res1@tmXBLabelStride = 4 ; else ; res1@trXMinF = iDOY ; res1@trXMaxF = iDOY+1 ; res1@tmXBValues = fspan(177,213,1+(213-177)*8) ; res1@tmXBLabels = (/"6/26","","6","","12","","18","","6/27","","6","","12","","18","","6/28","","6","","12","","18","","6/29","","6","","12","","18","","6/30","","6","","12","","18","","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 utc_date = cd_calendar(time, 0) delete(time) month = tointeger(utc_date(:,1)) ; use sprinti day = tointeger(utc_date(:,2)) delete(utc_date) ; hour = tointeger(utc_date(:,3)) ; year = tointeger(utc_date(:,0)) ; Convert to integer for ; minute = tointeger(utc_date(:,4)) ; second = utc_date(:,5) ; if(hour() day_0 = day(0) month_0 = month(0) day_1 = day_0+1 delete(month) delete(day) ; res1@tmXBMode = "Explicit" ; res1@tmXBMinorValues = ispan(0,48,3) ; res1@tmXBValues = ispan(6,48,6) ; UTC ; res1@tmXBLabels = (/month_0+"/"+day_0,"6","12","18",month_0+"/"+day_1,"6","12","18"/); CST res1@tmXBFormat = "f" res1@tmYLFormat = "f" ; res1@trXAxisType = "IrregularAxis";"LinearAxis" res1@tmYLMinorOn = False 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@tmYRMajorOutwardLengthF = 0 res1@tmYLMajorOutwardLengthF = 0 res1@tmXBMajorOutwardLengthF = 0 res1@tmXBMinorOutwardLengthF = 0 res1@cnFillOn = True res1@lbOrientation = "vertical" res1@lbTitleString = "WSP~C~ m s~S~-1~N~" res1@lbTitlePosition = "Top" ; title position res1@lbTitleDirection = "Across" ; title direction ; res1@lbTopMarginF = -0.05 res1@lbTopMarginF = 0.05 res1@lbLeftMarginF = -0.2 res1@lbRightMarginF = 0.8 res1@lbTitleOffsetF = 0.02 res1@cnLinesOn = False ; turn off contour lines res1@gsnAddCyclic = False res1@cnLineLabelsOn = False ; turn off contour labels res1@cnLevelSelectionMode = "ManualLevels" if(species_infile(ispecies).eq."WSPD")then res1@lbTitleOn = True if(year.eq.2003) then res1@cnMaxLevelValF = 20 else res1@cnMaxLevelValF = 24 end if res1@cnMinLevelValF = 1 res1@cnLevelSpacingF =1 else res1@lbTitleOn = False res1@cnMaxLevelValF = 350 res1@cnMinLevelValF = 10 res1@cnLevelSpacingF =10 end if res1@lbLabelStride =4 res1@lgPerimOn = False res1@lbLabelFontHeightF=0.015 res1@lgLabelFontHeightF = 0.0205 res1@lbTitleFontHeightF= .015 res1@tiXAxisFontHeightF= 0.0198 res1@tiYAxisFontHeightF= 0.0198 res1@tmXBLabelFontHeightF = 0.0198 res1@tmYLLabelFontHeightF = 0.0198 res1@cnFillMode = "RasterFill" index_missing=ind(ismissing(height)) delete(height) res1@gsnAddCyclic = False res1@gsnRightString = "Halo Lidar" ; print(index_missing(0)) printVarSummary(species_value) plot = gsn_csm_contour(wks, species_value(:,:), res1 ) ; draw second plot delete(index_missing) delete(species_value) ; system("convert -density 300 -trim -resize 800x800 "+figurename+".eps /nsftor/xhu/public_html/ARM/915rwpwindconC1/"+folder+"/"+figurename+".png") system("convert ARMlidartime_height_WSPD_neightCenter_0.eps ARMlidartime_height_WSPD_neightCenter_0.png") ; system("./convert.csh") ; system("ncl test_convert.ncl") ; system("rm *.eps ") ; system("mv "+figurename+".eps /nsftor/xhu/public_html/ARM/915rwpwindconC1/"+folder+"/.") print("finish "+ figurename) end do ; end do ; end do ; ispecies delete(ARMfiles) end do ; imonth end