load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" ;load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" ;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl" load "/home/carl/nclscripts/lib/tTick.ncl" load "/home/carl/nclscripts/lib/date2gem.ncl" load "/home/carl/nclscripts/lib/kf_filter.ncl" load "/home/carl/nclscripts/lib/ut_string.ncl" ;>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ;>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> begin print( "Making 260-day filtered VP plots..." ) plotType = "png" plotDir = "/ct13/ventrice/real_time/timeLon/" obsPerDay = 4 ; read the GFS 200-hPa winds print( "Reading 200 hPa Wind Data..." ) fileName = "/ct13/ventrice/real_time/gfs200Long.nc" inFile = addfile( fileName, "r" ) u = inFile->u200(:,:,:) v = inFile->v200(:,:,:) time = inFile->time sf = u vp = u uv2sfvpf(u,v,sf,vp) delete(sf) total = vp(:,:,:) print( "Missing zonal wind data: " + num(ismissing(total)) ) utc_date = ut_calendar( total&time, -5 ) ddd = day_of_year( utc_date(:,0), utc_date(:,1), utc_date(:,2) ) yyyyddd = ( 1000 * utc_date(:,0) ) + ddd print( "Calcuating the zonal wind anomalies..." ) print( "Reading the 200 hPa VP Climo..." ) fileName = "/ct13/ventrice/data/interim/raw/vp.200.climo.nc" inFile = addfile( fileName, "r" ) inVPclimo = lonFlip(inFile->climatology(:,:,:)) vpClimo = linint2_Wrap( inVPclimo&lon, inVPclimo&lat, inVPclimo, \\ True, v&lon, v&lat, 0 ) anom = calcDayAnomTLL( total, yyyyddd, vpClimo ) total = total * 10^-6 total = lonFlip(total) anom = anom * 10^-6 anom = lonFlip(anom) ;Start Padding with zeros ; ++++++++++++++++++++++++++ ; reduce some size And save some time on filtering ; ++++++++++++++++++++++++++ ; this saves 2/3's of the filtering computation time, ; which wasn't being used anyway. totalr = total(:,{-30:40},:) anomr = anom(:,{-30:40},:) delete(total) delete(anom) ; rename back to orig to work with rest of script total = totalr anom = anomr delete(totalr) delete(anomr) days = 120 days_of_zero = days * 2 * 4 ;change first number timeTot = dimsizes(time) - 1 stFill = days * 4 enFill = timeTot + stFill print(timeTot) print(stFill) print(enFill) total_dims = dimsizes(total) total_pad = new( (/ total_dims(0)+(days_of_zero) , total_dims(1), total_dims(2)/), typeof(total) ) anom_pad = total_pad total_pad = 0 total_pad( stFill:enFill,:,: ) = total ;480 because 360 + 120 for 30 after before anom_pad = 0 anom_pad( stFill:enFill,:,: ) = anom time_pad = fspan(total_pad&time(stFill)-days,total_pad&time(stFill),stFill+1) total_pad&time(:stFill) = time_pad time_pad = fspan(total_pad&time(enFill),total_pad&time(enFill)+days,stFill+1) total_pad&time(enFill:) = time_pad anom_pad&time = total_pad&time ;---------For filter------------- mjo_total = total_pad mjo_total@long_name = "Madden-Julian Oscillation" mjo_total@filter = "Wheeler & Kiladis (1999)" mjo_total@wavenumber = (/ 1, 5 /) mjo_total@period = (/ 30, 96 /) mjo_anom = anom_pad mjo_anom@long_name = "Madden-Julian Oscillation" mjo_anom@filter = "Wheeler & Kiladis (1999)" mjo_anom@wavenumber = (/ 1, 5 /) mjo_anom@period = (/ 30, 96 /) kelvin_total = total_pad kelvin_total@long_name = "Kelvin waves" kelvin_total@filter = "Wheeler & Kiladis (1999)" kelvin_total@depth = (/ 8, 90 /) kelvin_total@wavenumber = (/ 1, 14 /) kelvin_total@period = (/ 2.5, 30 /) kelvin_anom = anom_pad kelvin_anom@long_name = "Kelvin waves" kelvin_anom@filter = "Wheeler & Kiladis (1999)" kelvin_anom@depth = (/ 8, 90 /) kelvin_anom@wavenumber = (/ 1, 14 /) kelvin_anom@period = (/ 2.5, 30 /) er_total = total_pad er_total@long_name = "Equatorial Rossby waves" er_total@filter = "Wheeler & Kiladis (1999)" er_total@depth = (/ 8, 90 /) er_total@wavenumber = (/ -10, -1 /) er_total@period = (/ 5, 48 /) er_anom = anom_pad er_anom@long_name = "Equatorial Rossby waves" er_anom@filter = "Wheeler & Kiladis (1999)" er_anom@depth = (/ 8, 90 /) er_anom@wavenumber = (/ -10, -1 /) er_anom@period = (/ 5, 48 /) td_anom = anom_pad td_anom@long_name = "TD-type disturbances" td_anom@filter = "Roundy and Frank (2004), extended to wavenumber 20 (similar to Kiladis et al. 2006) and to 10 day waves for over Africa" td_anom@wavenumber = (/ -20, -6 /) td_anom@period = (/ 2, 10 /) mis = -999 mis@_FillValue = -999 print( "Filtering..." ) do y = 0, ( dimsizes(total_pad&lat) - 1 ) print( y + " " + ( dimsizes(total_pad&lat) - 1 ) ) kelvin_total(time|:,lat|y,lon|:) = (/ kf_filter( total_pad(time|:,lat|y,lon|:), \\ obsPerDay, 2.5, 30, 1, 14, 5, 90, "Kelvin" ) /) kelvin_anom(time|:,lat|y,lon|:) = (/ kf_filter( anom_pad(time|:,lat|y,lon|:), \\ obsPerDay, 2.5, 30, 1, 14, 5, 90, "Kelvin" ) /) er_total(time|:,lat|y,lon|:) = (/ kf_filter( total_pad(time|:,lat|y,lon|:), \\ obsPerDay, 05, 48, -10, -1, 8, 90, "ER" ) /) er_anom(time|:,lat|y,lon|:) = (/ kf_filter( anom_pad(time|:,lat|y,lon|:), \\ obsPerDay, 05, 48, -10, -1, 8, 90, "ER" ) /) mjo_total(time|:,lat|y,lon|:) = (/ kf_filter( total_pad(time|:,lat|y,lon|:), \\ obsPerDay, 30, 96, 1, 5, mis, mis, "None" ) /) mjo_anom(time|:,lat|y,lon|:) = (/ kf_filter( anom_pad(time|:,lat|y,lon|:), \\ obsPerDay, 30, 96, 1, 5, mis, mis, "None" ) /) td_anom(time|:,lat|y,lon|:) = (/ kf_filter( anom_pad(time|:,lat|y,lon|:), \\ obsPerDay, 2, 10, -20, -6, 90, mis, "mrg" ) /) end do print("done filtering") diro = "/free4/ventrice/" filo = "VP200_filtered_RT.nc" system("/bin/rm -f "+diro+filo) ; remove any pre-existing file ncdf = addfile(diro+filo,"c") ; open output netCDF file ; make time an UNLIMITED dimension filedimdef(ncdf,"time",-1,True) ; recommended for most applications ; output variables directly ncdf->kelvin_anom = kelvin_anom ncdf->er_anom = er_anom ncdf->mjo_anom = mjo_anom ncdf->td_anom = td_anom ;------------------ Read in STDs ------------------------------------------- kelvin_std_in = addfile("/free3/ventrice/data/interim/vp200/kelvin_stnd.nc","r") kelvin_stnd = kelvin_std_in->kelvin_stnd({-30:40},:) mjo_std_in = addfile("/free3/ventrice/data/interim/vp200/mjo_stnd.nc","r") mjo_stnd = mjo_std_in->mjo_stnd({-30:40},:) er_std_in = addfile("/free3/ventrice/data/interim/vp200/er_stnd.nc","r") er_stnd = er_std_in->er_stnd({-30:40},:) dim_time = dimsizes(kelvin_anom&time) lon = fspan(-180,180, 240) lat = fspan(-30,40, 47) kelvin_anoms = area_hi2lores_Wrap (kelvin_anom&lon,kelvin_anom&lat, kelvin_anom , True, 1, lon, lat, False) mjo_anoms = area_hi2lores_Wrap (mjo_anom&lon,mjo_anom&lat, mjo_anom , True, 1, lon, lat, False) er_anoms = area_hi2lores_Wrap (er_anom&lon,er_anom&lat, er_anom , True, 1, lon, lat, False) kelvin_totals = area_hi2lores_Wrap (kelvin_total&lon,kelvin_total&lat, kelvin_total , True, 1, lon, lat, False) mjo_totals = area_hi2lores_Wrap (mjo_total&lon,mjo_total&lat, mjo_total , True, 1, lon, lat, False) er_totals = area_hi2lores_Wrap (er_total&lon,er_total&lat, er_total , True, 1, lon, lat, False) print("Removing STD for Kelvin") do i = 0, dim_time-1 kelvin_anoms(i,:,:) = kelvin_anoms(i,:,:) / kelvin_stnd end do print("Removing STD for MJO") do i = 0, dim_time-1 mjo_anoms(i,:,:) = mjo_anoms(i,:,:) / mjo_stnd end do print("Removing STD for ERW") do i = 0, dim_time-1 er_anoms(i,:,:) = er_anoms(i,:,:) / er_stnd end do print("Removing STD for Kelvin") do i = 0, dim_time-1 kelvin_totals(i,:,:) = kelvin_totals(i,:,:) / kelvin_stnd end do print("Removing STD for MJO") do i = 0, dim_time-1 mjo_totals(i,:,:) = mjo_totals(i,:,:) / mjo_stnd end do print("Removing STD for ERW") do i = 0, dim_time-1 er_totals(i,:,:) = er_totals(i,:,:) / er_stnd end do res = True res@cnFillOn = True res@cnMonoFillColor = False res@cnLineLabelsOn = False res@cnInfoLabelOn = False res@cnLinesOn = False res@gsnSpreadColorStart = 2 res@gsnSpreadColorEnd = 21 res@gsnSpreadColors = True res@cnLevelSelectionMode = "ManualLevels" res@cnMinLevelValF = -15 ; set min contour level res@cnMaxLevelValF = 15 ; set max contour level res@cnLevelSpacingF = 1.5 ; set contour spacing res@cnMissingValFillColor = "gray75" res@cnMissingValFillPattern = 0 res@gsnRightString = "" ;res@gsnAddCyclic = True res@cnFillDrawOrder = "PreDraw" res@trYReverse = True res@gsnDraw = False res@gsnFrame = False res@lbTitleString = " [10~S~6~N~ m~S~2~N~s~S~-1~N~]" res@lbTitlePosition = "Right" res@lbTitleDirection = "Across" res@lbTitleFontHeightF = 0.02 res@lbBoxLinesOn = True ; res@lbTopMarginF = 0.5 ; res@pmLabelBarWidthF = 0.4 ;Make pretty month labels resTick = True resTick@Format = "%c %D" resTick@TickAxis = "YL" resTick@MajorStride = 80 ;30 resTick@Minorstride = 40 ;10 tTick( total&time, res, resTick ) anomRes = res ;anomRes@cnLevels = ispan( -8, 8, 1 ) * 1.5 totalRes = res ; totalRes@cnLevels = ispan( -8, 8, 1 ) * 1.5 ;Resources for KW Contours kwRes = res kwRes@gsnRightString = " " kwRes@gsnLeftString = " " kwRes@cnFillOn = False kwRes@cnLinesOn = True kwRes@cnLevelSelectionMode = "ExplicitLevels" kwRes@cnLevels = (/-4,-2,-1,1,2,4/) kwRes@cnLineThicknessF = 1 kwRes@gsnContourPosLineDashPattern = 0 kwRes@gsnContourNegLineDashPattern = 3 kwRes@cnLineColor = "black" ;Resources for MJO Contours mjoRes = res mjoRes@gsnRightString = " " mjoRes@gsnLeftString = " " mjoRes@cnFillOn = False mjoRes@cnLinesOn = True mjoRes@cnLevelSelectionMode = "ExplicitLevels" mjoRes@cnLevels = (/ -5,-4.5,-4,-3.5,-3,-2.5,-2,-1.5,-1,-0.5,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5/) mjoRes@cnLineThicknessF = 1.5 mjoRes@gsnContourPosLineDashPattern = 0 mjoRes@gsnContourNegLineDashPattern = 3 mjoRes@cnLineColor = "red" ;Resources for ER Contours erRes = res erRes@gsnRightString = " " erRes@gsnLeftString = " " erRes@cnFillOn = False erRes@cnLinesOn = True erRes@cnLevelSelectionMode = "ExplicitLevels" erRes@cnLevels = (/ -4,-2,2,4/) erRes@cnLineThicknessF = 1.5 erRes@gsnContourPosLineDashPattern = 0 erRes@gsnContourNegLineDashPattern = 3 erRes@cnLineColor = "darkgreen" if( ( plotType.eq."png" ).or.( plotType.eq."gif" ) ) then plotTypeLocal = "eps" else plotTypeLocal = plotType end if minLat = (/ -30, -20, -10, -5, 0, 10, 20 /) maxLat = (/ -20, -10, 0, 5, 10, 20, 30 /) do i = 0, dimsizes( minLat)-1 latLabel = abs( minLat(i) ) + where( minLat(i).lt.0, "S", "N" ) + "-" \\ + abs( maxLat(i) ) + where( maxLat(i).le.0, "S", "N" ) plotName = "vp.filtLong.total.260." + latLabel print( (/ plotName /) ) plotName = plotDir + plotName totalRes@gsnLeftString = latLabel + ": 200-hPa velocity potential" xtData = dim_avg_n_Wrap( total(:,{minLat(i):maxLat(i)},:),1) wrf_smooth_2d(xtData(:,:),2) kwData = dim_avg_n_Wrap( kelvin_totals(:,{minLat(i):maxLat(i)},:),1) mjoData = dim_avg_n_Wrap( mjo_totals( :,{minLat(i):maxLat(i)},:),1) ; erData = dim_avg_n_Wrap( er_totals(:,{minLat(i):maxLat(i)},:),1) wks = gsn_open_wks( plotTypeLocal, plotName ) gsn_merge_colormaps(wks,"vp","default") plot = gsn_csm_hov( wks, xtData, totalRes ) kwPlot = gsn_csm_hov( wks, kwData, kwRes ) mjoPlot = gsn_csm_hov( wks, mjoData, mjoRes ) ; erPlot = gsn_csm_hov( wks, erData, erRes ) overlay(plot, mjoPlot) ;overlay(plot, erPlot) overlay(plot, kwPlot) draw(plot) frame(wks) if( ( plotType.eq."png" ).or.( plotType.eq."gif" ) ) then system( "convert -trim -density 200 " \\ + plotName + ".eps " + plotName + "." + plotType ) system( "rm -f " + plotName + ".eps" ) end if plotName = "vp.filtLong.anom.260." + latLabel print( (/ plotName /) ) plotName = plotDir + plotName anomRes@gsnLeftString = latLabel + ": 200-hPa velocity potential anomalies" xtData = dim_avg_n_Wrap( anom(:,{minLat(i):maxLat(i)},:),1) wrf_smooth_2d(xtData(:,:),2) kwData = dim_avg_n_Wrap( kelvin_anoms(:,{minLat(i):maxLat(i)},:),1) mjoData = dim_avg_n_Wrap( mjo_anoms(:,{minLat(i):maxLat(i)},:),1) ; erData = dim_avg_n_Wrap( er_anoms(:,{minLat(i):maxLat(i)},:),1) wks = gsn_open_wks( plotTypeLocal, plotName ) gsn_merge_colormaps(wks,"vp","default") plot = gsn_csm_hov( wks, xtData, anomRes ) kwPlot = gsn_csm_hov( wks, kwData, kwRes ) mjoPlot = gsn_csm_hov( wks, mjoData, mjoRes ) ; erPlot = gsn_csm_hov( wks, erData, erRes ) overlay(plot, mjoPlot) ; overlay(plot, erPlot) overlay(plot, kwPlot) draw(plot) frame(wks) if( ( plotType.eq."png" ).or.( plotType.eq."gif" ) ) then system( "convert -trim -density 200 " \\ + plotName + ".eps " + plotName + "." + plotType ) system( "rm -f " + plotName + ".eps" ) end if end do print( "" ) print( "Thank you, come again." ) end ;<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ;<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<