;*********************************************************** ; Combined EOFs ;*********************************************************** 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/csm/contributed.ncl" begin neof = 2 latS = -15 latN = 15 plotType = "png" plotName = "kelvin_vp200_clivar" plotDir = "/ct13/ventrice/real_time/timeLon/" ; plot directory diri = "/free4/ventrice/" filu200 = "VP200_filtered_RT.nc" fileName = "/free4/ventrice/VP200_filtered_RT.nc" inFile = addfile( fileName, "r" ) kelvin = lonFlip(inFile->kelvin_anom(442:1550,:,:) ) print( "Missing zonal wind data: " + num(ismissing(kelvin)) ) utc_date = ut_calendar( kelvin&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..." ) lon = fspan(-180,180, 240) lat = fspan(-90,90, 121) KELVIN = area_hi2lores_Wrap (kelvin&lon,kelvin&lat, kelvin , True, 1, lon, lat, False) time = kelvin&time KELVIN!0 = "time" KELVIN!1 = "lat" KELVIN!2 = "lon" ;*********************************************************** ; Read anomalies ;*********************************************************** vp_1 = KELVIN(lon|:,time|:, lat|:) VP = dim_avg_Wrap(vp_1(:,:,{latS:latN})) dimw = dimsizes( KELVIN ) ntim = dimw(0) nlat = dimw(1) mlon = dimw(2) ;************************************************ ; Compute the temporal variance ;************************************************ var_vp = dim_variance_Wrap( VP) ; (lon) ;************************************************ ; Compute the zonal mean of the temporal variance ;************************************************ zavg_var_vp = dim_avg_Wrap( var_vp ) ;************************************************ ; Normalize by sqrt(avg_var*) ;************************************************ VP = VP/sqrt(zavg_var_vp ) ; (lon,time) ;************************************************ ; Combine the normalized data into one variable ;************************************************ printVarSummary(VP) cdata = new ( (/mlon,ntim/), typeof(VP), getFillValue(VP)) printVarSummary(cdata) do ml=0,mlon-1 cdata(ml ,:) = (/ VP(ml,:) /) end do ;************************************************ ; Compute Combined EOF ;************************************************ ; eof_cdata = eofunc_Wrap(cdata , neof, False) ; (neof,3*mlon) fileName = "/free4/ventrice/kelvin/filterVP/eof_cdata.nc" inFile = addfile( fileName, "r" ) eof_cdata = inFile->eof_cdata eof_ts_cdata = eofunc_ts_Wrap(cdata,eof_cdata,False) ; (neof,time) print("==============") printVarSummary(eof_cdata) printMinMax(eof_cdata, True) print("==============") printVarSummary(eof_ts_cdata) printMinMax(eof_ts_cdata, True) ;************************************************ ; For clarity, explicitly extract each variable ;************************************************ nvar = 1 ; "vp" ceof = new( (/nvar,neof,mlon/), typeof(cdata), getFillValue(cdata)) do n=0,neof-1 ceof(0,n,:) = eof_cdata(n,0:mlon-1) ; vp end do ceof_ts = new( (/nvar,neof,ntim/), typeof(cdata), getFillValue(cdata)) ceof_ts(0,:,:) = eofunc_ts_Wrap( VP,ceof(0,:,:),False) ; (neof,time) ;************************************************ ; Compute cross correlation of each variable's EOF time series at zero-lag ;************************************************ r_vp_vp = escorc(ceof_ts(0,:,:) , ceof_ts(0,:,:)) print("==============") do n=0,neof-1 print("neof="+n \ +" r_vp_vp=" +sprintf("%4.3f",r_vp_vp(n)) ) end do print("==============") ;************************************************ ; Compute cross correlation of the multivariate EOF; EOF 1 vs EOF 2 ;************************************************ mxlag = 25 rlag_01 = esccr(eof_ts_cdata(0,:),eof_ts_cdata(0,:), mxlag) ; (N,mxlag+1) rlag_10 = esccr(eof_ts_cdata(0,:),eof_ts_cdata(0,:), mxlag) ; (N,mxlag+1) ccr_12 = new ( (/2*mxlag+1/), float) ccr_12(mxlag:) = rlag_10(0:mxlag) ccr_12(0:mxlag) = rlag_01(::-1) ; reverse order ;************************************************ ; Normalize the multivariate EOF 1&2 component time series ; Compute (PC1^2+PC2^2): values > 1 indicate "strong" periods ;************************************************ eof_ts_cdata(0,:) = eof_ts_cdata(0,:)/stddev(eof_ts_cdata(0,:)) eof_ts_cdata(1,:) = eof_ts_cdata(1,:)/stddev(eof_ts_cdata(1,:)) ; eof_ts_cdata(1,:) = -1. * eof_ts_cdata(1,:) ; eof_ts_cdata(0,:) = -1. * eof_ts_cdata(0,:) kelvin_ts_index = eof_ts_cdata(0,:)^2 + eof_ts_cdata(1,:)^2 kelvin_ts_index_smt = runave(kelvin_ts_index, 91, 0) ; 91-day running mean nGood = num(.not.ismissing(kelvin_ts_index)) ; # non-missing nStrong = num(kelvin_ts_index .ge. 1.0) print("nGood="+nGood+" nStrong="+nStrong+" nOther="+(nGood-nStrong)) ;************************************************ ; Write PC results to netCDF for use in another example. ;************************************************ kelvin_ts_index!0 = "time" kelvin_ts_index&time = time kelvin_ts_index@long_name = "KELVIN PC INDEX" kelvin_ts_index@info = "(PC1^2 + PC2^2)" PC1 = eof_ts_cdata(0,:) PC1!0= "time" PC1&time = time PC1@long_name = "PC1" PC1@info = "PC1/stddev(PC1)" PC2 = eof_ts_cdata(1,:) PC2!0= "time" PC2&time = time PC2@long_name = "PC2" PC2@info = "PC2/stddev(PC2)" diro = "/free4/ventrice/kelvin/filterVP/" filo = "KELVIN_VP200_PC_INDEX.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_INDEX = kelvin_ts_index ncdf->PC1 = PC1 ncdf->PC2 = PC2 ;------------------------------------------------------------ ; PLOTS ;------------------------------------------------------------ yyyymmdd = ut_calendar(time, -3) yrfrac = yyyymmdd_to_yyyyfrac(yyyymmdd, 0.0) delete(yrfrac@long_name) delete(lon@long_name) day = ispan(-mxlag, mxlag, 1) ;day@long_name = "lag (day)" if (plotType.eq."png") then pltTypeLocal = "eps" else pltTypeLocal = plotType end if pltPath = "/ct13/ventrice/real_time/"+plotName wks = gsn_open_wks(pltTypeLocal,pltPath) plot = new(3,graphic) ;************************************************ ; Multivariate EOF plots ;************************************************ rts = True rts@gsnDraw = False ; don't draw yet rts@gsnFrame = False ; don't advance frame yet rts@gsnScale = True ; force text scaling rts@vpHeightF = 0.40 ; Changes the aspect ratio rts@vpWidthF = 0.85 rts@vpXF = 0.10 ; change start locations rts@vpYF = 0.75 ; the plot rts@xyLineThicknesses = (/2, 2, 2/) rts@gsnYRefLine = 0. ; reference line rts@pmLegendDisplayMode = "Always" ; turn on legend rts@pmLegendSide = "Top" ; Change location of rts@pmLegendParallelPosF = 0.86 ; move units right rts@pmLegendOrthogonalPosF = -0.50 ; move units down rts@pmLegendWidthF = 0.15 ; Change width and rts@pmLegendHeightF = 0.15 ; height of legend. rts@lgLabelFontHeightF = 0.0175 rtsP = True ; modify the panel plot rtsP@gsnMaximize = True ; large format rtsP@txString = "Multivariate EOF: 15S-15N - Kelvin filtered VP200" rts@xyExplicitLegendLabels = (/"VP"/) do n=0,neof-1 rts@gsnLeftString = "EOF "+(n+1) rts@gsnRightString = sprintf("%3.1f",ceof@pcvar(n)) +"%" plot(n) = gsn_csm_xy (wks,lon,ceof(:,n,:),rts) end do gsn_panel(wks,plot(0:1),(/2,1/),rtsP) ; now draw as one plot ;************************************************ ; ? png ? ;************************************************ if( ( plotType.eq."png" ).or.( plotType.eq."gif" ) ) then system( "convert -trim +repage -density 144 " \\ + pltPath + ".eps " + pltPath + "." + plotType ) system( "rm -f " + pltPath + ".eps" ) end if end