load "~/loadFiles.ncl" ; load "/nfs/lb13/ppapin/scripts/ncl/circ/radavg/radavg.ncl" ; load "/home11/grad/2010/abrammer/NCL/functions/uv2cvF_Wrap.ncl" diri = ".." diro = "./" ; diro = "/nfs/ct12/yc689758/ATM741/AEW-TC/cross-sec/WITH_VORTEX_27km" latRange=(/0,20/) fin = "TEST_wrfout_d03_2004-08-11_18:00:00.nc" inf = addfile(diri+"/"+fin,"r") ;---Get input file to get lat,lon dimensions lat = inf->XLAT(0,:,:) lon = inf->XLONG(0,:,:) ntimes = tointeger(systemfunc("ls -1R "+diri+"/TEST* | wc -l")) ;---What times and how many time steps are in the data set? year = "2004" month = "08" day = "10" hour = "18" type = "png" do it = 0,ntimes-1,1 ;---Time loop every 6 hours ; do it = 4,ntimes-1,1 ;---Time loop every 6 hours print("Working on: "+year+"-"+month+"-"+day+"_"+hour) a := addfile(diri+"/TEST_wrfout_d03_"+year+"-"+month+"-"+day+"_"+hour+":00:00.nc","r") print("reading in files") p_base := wrf_user_getvar(a,"PB",0) p_prime := wrf_user_getvar(a,"P",0) p := p_base + p_prime u := wrf_user_getvar(a,"ua",0) ;10=925mb, 16=600mb v := wrf_user_getvar(a,"va",0) ;10=925mb, 16=600mb ; avo:= wrf_user_getvar(a,"avo",0) ;10=925mb, 16=600mb ; pv := wrf_user_getvar(a,"pvo",0) ;10=925mb, 16=600mb rh := wrf_user_getvar(a,"rh",0) ; relative humidity qv := wrf_user_getvar(a,"QVAPOR",0) ; relative humidity dumlevs = ispan(100,1000,50) levs = dumlevs(::-1) ndim=dimsizes(u) var_u:=new((/dimsizes(levs),ndim(1),ndim(2)/),typeof(u)) var_u!0 = "lev" var_u&lev = levs var_u&lev@units = "hPa" var_u!1 = "lat" var_u&lat = (/lat(:,0)/) var_u&lat@units="degrees_north" var_u&lat@description="latitude" var_u!2 = "lon" var_u&lon = (/lon(0,:)/) var_u&lon@units="degrees_east" var_u&lon@description="longitude" var_v := var_u var_vort := var_u ; var_pv:=var_u ; var_avo:=var_u var_rh:=var_u var_qv:=var_u print("interpolating...") var_u = (/wrf_user_intrp3d(u,p,"h",levs,0.,False)/) var_v = (/wrf_user_intrp3d(v,p,"h",levs,0.,False)/) var_vort = (/uv2vr_cfd(var_u,var_v,lat(:,0),lon(0,:),0)*10^5/) ; var_avo = (/wrf_user_intrp3d(avo,p,"h",levs,0.,False)*10/) var_rh = (/wrf_user_intrp3d(rh,p,"h",levs,0.,False)/) var_qv = (/wrf_user_intrp3d(qv,p,"h",levs,0.,False)*10^3/) ;-----PERT POTENTIAL TEMPERATURE---------------------------------------------------------- ; lonRange=(/0,-5,-10,-15,-20,-25,-30/) lonRange=(/-10/) do i=0, dimsizes(lonRange)-1 print("vort") res = True res@gsnDraw = False res@gsnFrame = False res@cnFillOn = True res@cnLinesOn = False res@lbLabelFontHeightF = 0.025 res@tiXAxisFontHeightF = 0.020 res@tiYAxisFontHeightF = 0.020 res@gsnStringFontHeightF = 0.020 res@tmXBLabelFontHeightF = 0.016 res@tmYLLabelFontHeightF = 0.016 res@gsnRightString = hour+"00"+" UTC "+day+" Aug 2004" res@lbOrientation = "Vertical" res@pmTickMarkDisplayMode = "Always" ; turn on tickmarks res@tfDoNDCOverlay = True ; do not transform data res@gsnAddCyclic = False colorbar = "temp_diff_18lev" res@cnLevelSelectionMode = "ManualLevels" ; set manual contour l$ res@cnMinLevelValF = -10. ;-9 ; set min contou$ res@cnMaxLevelValF = 10. ;9 ; set max contour le$ res@cnLevelSpacingF = 1. ;1 ; set contour spacing res@tmYROn = False ; if(i .le. 5) res@gsnLeftString = "Vorticity ("+abs(lonRange(i))+"~S~o~N~W, 10~S~-5~N~ s~S~-1~N~)" ; units: 10-1 PVU res@tmYRMode = "Automatic" res@trYReverse = True system("mkdir -p "+ diro+"vort/") name = diro+"vort/vort_"+lon_convert(lonRange(i))+"_"+it wks = gsn_open_wks(type,name) gsn_define_colormap(wks,colorbar) plot = gsn_csm_pres_hgt(wks,var_vort({:200},{latRange(0):latRange(1)},{lonRange(i)}),res) draw(plot) frame(wks) ; end if ; system("convert -background white -flatten -trim +repage -density 100x100 -size 1000x1500 "+name+".png -trim +repage "+name+".png") ; ;-----RELATIVE HUMIDITY------------------------------------------------------------------- print("RH") res@cnLevelSelectionMode = "ManualLevels" ; set manual contour l$ res@cnMinLevelValF = 0. ;-9 ; set min contou$ res@cnMaxLevelValF = 100. ;9 ; set max contour le$ res@cnLevelSpacingF = 5 ;1 ; set contour spacing res@gsnLeftString = "Relative Humidity ("+abs(lonRange(i))+"~S~o~N~W)" system("mkdir -p "+ diro+"rh/") name = diro+"rh/rh_"+lon_convert(lonRange(i))+"_"+it wks9 = gsn_open_wks(type,name) gsn_define_colormap(wks9,"GMT_drywet") plot = gsn_csm_pres_hgt(wks9,var_rh({:200},{latRange(0):latRange(1)},{lonRange(i)}),res) draw(plot) frame(wks9) ; system("convert -background white -flatten -trim +repage -density 100x100 -size 1000x1500 "+name+".png -trim +repage "+name+".png") ; ;-----U WIND------------------------------------------------------------------------------ print("u") res@cnLevelSelectionMode = "ManualLevels" ; set manual contour l$ res@cnMinLevelValF = -12. ;-9 ; set min contou$ res@cnMaxLevelValF = 12. ;9 ; set max contour le$ res@cnLevelSpacingF = 1 ;1 ; set contour spacing res@gsnLeftString = "u ("+abs(lonRange(i))+"~S~o~N~W)" system("mkdir -p "+ diro+"u/") name = diro+"u/u_"+lon_convert(lonRange(i))+"_"+it wks10 = gsn_open_wks(type,name) colorbar = "temp_diff_18lev" gsn_define_colormap(wks10,colorbar) contour = gsn_csm_pres_hgt(wks10,var_u({:200},{latRange(0):latRange(1)},{lonRange(i)}),res) draw(contour) frame(wks10) ; system("convert -background white -flatten -trim +repage -density 100x100 -size 1000x1500 "+name+".png -trim +repage "+name+".png") ; ; ;-----V WIND------------------------------------------------------------------------------ print("v") res@cnLevelSelectionMode = "ManualLevels" ; set manual contour l$ res@cnMinLevelValF = -12. ;-9 ; set min contou$ res@cnMaxLevelValF = 12. ;9 ; set max contour le$ res@cnLevelSpacingF = 1 ;1 ; set contour spacing res@gsnLeftString = "v ("+abs(lonRange(i))+"~S~o~N~W)" system("mkdir -p "+ diro+"v/") name = diro+"v/v_"+lon_convert(lonRange(i))+"_"+it wks11 = gsn_open_wks(type,name) gsn_define_colormap(wks11,colorbar) contour = gsn_csm_pres_hgt(wks11,var_v({:200},{latRange(0):latRange(1)},{lonRange(i)}),res) draw(contour) frame(wks11) ; system("convert -background white -flatten -trim +repage -density 100x100 -size 1000x1500 "+name+".png -trim +repage "+name+".png") ; ;-----WATER VAPOR------------------------------------------------------------------------- print("qv") res@cnLevelSelectionMode = "ManualLevels" ; set manual contour l$ res@cnMinLevelValF = 1. ;-9 ; set min contou$ res@cnMaxLevelValF = 16. ;9 ; set max contour le$ res@cnLevelSpacingF = 1 ;1 ; set contour spacing res@gsnLeftString = "Mixing Ratio ("+abs(lonRange(i))+"~S~o~N~W)" system("mkdir -p "+ diro+"qv/") name = diro+"qv/qv_"+lon_convert(lonRange(i))+"_"+it ;units: g/kg wks15 = gsn_open_wks(type,name) colorbar = "GMT_drywet" gsn_define_colormap(wks15,colorbar) contour = gsn_csm_pres_hgt(wks15,var_qv({:200},{latRange(0):latRange(1)},{lonRange(i)}),res) draw(contour) frame(wks15) ; system("convert -background white -flatten -trim +repage -density 100x100 -size 1000x1500 "+name+".png -trim +repage "+name+".png") ;-----RELATIVE VORTICITY------------------------------------------------------------------ ; vort_600 = vort_600*10^5 ; vort_600@long_name = "Relative Vorticity (600 hPa)" ; res@cnLevelSelectionMode = "ManualLevels" ; set manual contour l$ ; res@cnMinLevelValF = -10. ;-9 ; set min contou$ ; res@cnMaxLevelValF = 10. ;9 ; set max contour le$ ; res@cnLevelSpacingF = 1 ;1 ; set contour spacing ; name = "vort_600_"+it ; wks12 = gsn_open_wks(type,name) ; gsn_define_colormap(wks12,colorbar) ; contour = gsn_csm_contour_map(wks12,vort_600(:,:),res) ; draw(contour) ; frame(wks12) ; system("convert -background white -flatten -trim +repage -density 100x100 -size 1000x1500 "+name+".png -trim +repage "+name+".png") ; ;-----CURVATURE VORTICITY----------------------------------------------------------------- ; curv_vort_600@long_name = "Curvature Vorticity (600 hPa)" ; res@cnLevelSelectionMode = "ManualLevels" ; set manual contour l$ ; res@cnMinLevelValF = -10 ;-9 ; set min contou$ ; res@cnMaxLevelValF = 10 ;9 ; set max contour le$ ; res@cnLevelSpacingF = 1 ;1 ; set contour spacing ; name = "curv_vort_600_"+it ; wks13 = gsn_open_wks(type,name) ; gsn_define_colormap(wks13,colorbar) ; contour = gsn_csm_contour_map(wks13,curv_vort_600(:,:),res) ; draw(contour) ; frame(wks13) ; system("convert -background white -flatten -trim +repage -density 100x100 -size 1000x1500 "+name+".png -trim +repage "+name+".png") ; ; ;-----CIRCULATION------------------------------------------------------------------------- ; circ_600 = circ_600*10^5 ; circ_600@long_name = "Circulation (600 hPa)" ; res@cnLevelSelectionMode = "ManualLevels" ; set manual contour l$ ; res@cnMinLevelValF = -10. ;-9 ; set min contou$ ; res@cnMaxLevelValF = 10. ;9 ; set max contour le$ ; res@cnLevelSpacingF = 1 ;1 ; set contour spacing ; name = "circ_600_"+it ; wks14 = gsn_open_wks(type,name) ; gsn_define_colormap(wks14,colorbar) ; contour = gsn_csm_contour_map(wks14,circ_600(:,:),res) ; draw(contour) ; frame(wks14) ; system("convert -background white -flatten -trim +repage -density 100x100 -size 1000x1500 "+name+".png -trim +repage "+name+".png") delete(var_rh) delete(var_u) delete(var_v) delete(var_vort) delete(var_qv) end do if(hour.eq."00") then hour = "01" else if(hour.eq."01") then hour = "02" else if(hour.eq."02") then hour = "03" else if(hour.eq."03") then hour = "04" else if(hour.eq."04") then hour = "05" else if(hour.eq."05") then hour = "06" else if(hour.eq."06") then hour = "07" else if(hour.eq."07") then hour = "08" else if(hour.eq."08") then hour = "09" else if(hour.eq."09") then hour = "10" else if(hour.eq."10") then hour = "11" else if(hour.eq."11") then hour = "12" else if(hour.eq."12") then hour = "13" else if(hour.eq."13") then hour = "14" else if(hour.eq."14") then hour = "15" else if(hour.eq."15") then hour = "16" else if(hour.eq."16") then hour = "17" else if(hour.eq."17") then hour = "18" else if(hour.eq."18") then hour = "19" else if(hour.eq."19") then hour = "20" else if(hour.eq."20") then hour = "21" else if(hour.eq."21") then hour = "22" else if(hour.eq."22") then hour = "23" else if(hour.eq."23") then hour = "00" day_int = tointeger(day) if(day_int.le.8) then day = "0"+tostring(day_int+1) else day = tostring(day_int+1) end if end if end if end if end if end if end if end if end if end if end if end if end if end if end if end if end if end if end if end if end if end if end if end if end if end do ; END OF TIME LOOP ; res = True ; res@gsnDraw = False ; res@gsnFrame = False ; ; res@gsnLeftString = " " ; turn off strings ; ; res@gsnRightString = " " ; turn off strings ; res@cnFillOn = True ; ; res@cnLinesOn = False ; res@cnInfoLabelOn = False ; res@cnLevelSelectionMode = "ExplicitLevels" ; cmap = RGBtoCmap("../sunshine_diff_20lev.2.rgb") ; res@cnLevelSelectionMode = "ExplicitLevels" ; res@cnFillColors = cmap(2::,:) ; ; res@cnLevels = (/-4.,-2.5,-1.75,-1,-.5,0,.5,1.,1.25,1.5,1.75,2,2.5,3.,3.5,4,4.5/)*2 ; res@cnLabelBarEndStyle = "ExcludeOuterBoxes" ; res@cnLevels = (/-5.,-4,-3,-2,-1.,0,1,2,3,4,5,6,7,8,9,10/)/4 ; res@lbLabelAutoStride = False ; res@lbBoxSeparatorLinesOn = False ; res@lbLabelFontHeightF = 0.015 ; ;res@trYReverse = True ; ; print("ploting.....") ; wcStrtClmP = systemfunc("date") ; plot = gsn_csm_pres_hgt(wks,var_3d({:200},50,:),res) ; ; I've also tried to change var_3d from (lev,lat,lon) to (lon,lev,lat) and see if ; ; it will speed up but it didn't work. ; ; The time it takes is also variable dependent, ua is faster than pvo ; ; In addition, res@cnFillMode = "CellFill" is faster than RasterFill and the default... ; wallClockElapseTime(wcStrtClmP, "Processing", 0) ; draw(plot) ; frame(wks)