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" lonRange=(/-40,10/) 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("get vars") 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/) ; creating the 3D coordinates for var_3d ; lat1d=lat(1,0) ; lon1d=lon(0,:) ; var_vort = uv2vr_cfd(var_u,var_v,lat1d,lon1d,0) ; var_avo=var_vort ;-----CALCULATE CURVATURE VORT AND CIRCULATION-------------------------------------------- ; rad = 100. ;---radius of circulation average around vortex; example in km ; wrf_resolution = 27. ;---wrf resultion in km (used in circulation calculation) ; mv = pv@_FillValue ; var_circ := var_pv ;-----circulation same dimensions at vort ; var_curv := var_pv ;-----circulation same dimensions at vort ; var_circ = 0. ; var_curv = 0. ; box = toint((rad / wrf_resolution)+2) ; var_curv = uv2cvF_Wrap(var_u,var_v) ; radavg(vort_925(:,:),circ_925(:,:), rad, box,mv); seems to be 2D only ;-----PERT POTENTIAL TEMPERATURE---------------------------------------------------------- ; lonRange=(/0,-5,-10,-15,-20,-25,-30/) latRange=(/5/) do i=0, dimsizes(latRange)-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 latRange@degree=True res@gsnLeftString = "Vorticity ("+lat_convert(latRange(i))+")" delete(latRange@degree) ; units: 10-1 PVU res@tmYRMode = "Automatic" res@trYReverse = True system("mkdir -p "+ diro+"vort/") name = diro+"vort/vort_"+lat_convert(latRange(i))+"_"+it wks = gsn_open_wks(type,name) gsn_define_colormap(wks,colorbar) plot = gsn_csm_pres_hgt(wks,var_vort({:200},{latRange(i)},{lonRange(0):lonRange(1)}),res) draw(plot) frame(wks) print("rh") ; system("convert -background white -flatten -trim +repage -density 100x100 -size 1000x1500 "+name+".png -trim +repage "+name+".png") ; ;-----RELATIVE HUMIDITY------------------------------------------------------------------- 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 latRange@degree=True res@gsnLeftString = "Relative Humidity ("+lat_convert(latRange(i))+")" delete(latRange@degree) system("mkdir -p "+ diro+"rh/") name = diro+"rh/rh_"+lat_convert(latRange(i))+"_"+it wks9 = gsn_open_wks(type,name) gsn_define_colormap(wks9,"GMT_drywet") plot = gsn_csm_pres_hgt(wks9,var_rh({:200},{latRange(i)},{lonRange(0):lonRange(1)}),res) draw(plot) frame(wks9) print("u") ; system("convert -background white -flatten -trim +repage -density 100x100 -size 1000x1500 "+name+".png -trim +repage "+name+".png") ; ;-----U WIND------------------------------------------------------------------------------ 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 latRange@degree=True res@gsnLeftString = "u ("+lat_convert(latRange(i))+")" delete(latRange@degree) system("mkdir -p "+ diro+"u/") name = diro+"u/u_"+lat_convert(latRange(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(i)},{lonRange(0):lonRange(1)}),res) draw(contour) frame(wks10) print("v") ; system("convert -background white -flatten -trim +repage -density 100x100 -size 1000x1500 "+name+".png -trim +repage "+name+".png") ; ; ;-----V WIND------------------------------------------------------------------------------ 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 latRange@degree=True res@gsnLeftString = "v ("+lat_convert(latRange(i))+")" delete(latRange@degree) system("mkdir -p "+ diro+"v/") name = diro+"v/v_"+lat_convert(latRange(i))+"_"+it wks11 = gsn_open_wks(type,name) gsn_define_colormap(wks11,colorbar) contour = gsn_csm_pres_hgt(wks11,var_v({:200},{latRange(i)},{lonRange(0):lonRange(1)}),res) draw(contour) frame(wks11) print("qv") ; system("convert -background white -flatten -trim +repage -density 100x100 -size 1000x1500 "+name+".png -trim +repage "+name+".png") ; ;-----WATER VAPOR------------------------------------------------------------------------- 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 latRange@degree=True res@gsnLeftString = "Mixing Ratio ("+lat_convert(latRange(i))+")" delete(latRange@degree) system("mkdir -p "+ diro+"qv/") name = diro+"qv/qv_"+lat_convert(latRange(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(i)},{lonRange(0):lonRange(1)}),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") 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)