;*********************************************************** ; 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 = "kelvinclivar" plotDir = "/ct13/ventrice/real_time/timeLon/" ; plot directory diri = "/ct13/ventrice/real_time/" filu200 = "gfs200.nc" fileName = "/ct13/ventrice/real_time/gfs200.nc" inFile = addfile( fileName, "r" ) u = inFile->u200(:,{-30:30},:) v = inFile->v200(:,{-30:30},:) sf = u vp = u uv2sfvpf(u,v,sf,vp) delete(sf) ;Test Matt's real time velocity potential ; diric = "/ct13/janiga/rtmaps/data2/vp.200.cat.nc" ; inMfile = addfile( diric, "r") ; velPot = inMfile->vp(:,{-30:30},:) total = vp printVarSummary(total) 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 = "/free3/ventrice/data/interim/raw/vp.200.climo.nc" fileName = "/ct13/janiga/data/era_interim/annual_cycle/vp.200.ann_cyc.1989-2009.nc" inFile = addfile( fileName, "r" ) ;inVPclimo = lonFlip(inFile->climatology({ddd},{-30:30},:)) inVPclimo = lonFlip(inFile->vp({ddd},{-30:30},:)) vpClimo = linint2_Wrap( inVPclimo&lon, inVPclimo&lat, inVPclimo, \\ True, v&lon, v&lat, 0 ) printVarSummary(total) printVarSummary(yyyyddd) printVarSummary(vpClimo) VPanom = calcDayAnomTLL( total, yyyyddd, vpClimo ) print( "Calcuating the 850 hPa zonal wind anomalies..." ) fileName = "/ct13/ventrice/real_time/gfs850.nc" inFile = addfile( fileName, "r" ) u850total = inFile->u850(:,{-30:30},:) fileName = "/ct13/ventrice/real_time/wind850.climo.nc" inFile = addfile( fileName, "r" ) if( ismissing( inFile ) ) then print( "DANGER! DANGER! File could not be opened: " + fileName ) system( "mail -s 'Error opening file: " + fileName \\ + "' ventrice@atmos.albany.edu < error_email.txt " ) end if in850Climo = inFile->u850(:,{-30:30},:) climo_850 = linint2_Wrap( in850Climo&lon, in850Climo&lat, in850Climo, \\ True, u850total&lon, u850total&lat, 0 ) u850anom = calcDayAnomTLL( u850total, yyyyddd, climo_850 ) print( "Calcuating the 200 hPa zonal wind anomalies..." ) fileName = "/ct13/ventrice/real_time/gfs200.nc" inFile = addfile( fileName, "r" ) u200total = inFile->u200(:,{-30:30},:) fileName = "/ct13/ventrice/data/interim/raw/u.200.climo.nc" inFile = addfile( fileName, "r" ) if( ismissing( inFile ) ) then print( "DANGER! DANGER! File could not be opened: " + fileName ) system( "mail -s 'Error opening file: " + fileName \\ + "' ventrice@atmos.albany.edu < error_email.txt " ) end if in200Climo = lonFlip(inFile->climatology(:,{-30:30},:)) climo_200 = linint2_Wrap( in200Climo&lon, in200Climo&lat, in200Climo, \\ True, u200total&lon, u200total&lat, 0 ) u200anom = calcDayAnomTLL( u200total, yyyyddd, climo_200 ) pltType = "eps" pltName = "kelvinclivar" ; yrStrt+"_"+yrLast ;*********************************************************** ; Read anomalies ;*********************************************************** vp_1 = VPanom(lon|:,time|:, lat|:) VP = dim_avg_Wrap(vp_1(:,:,{latS:latN})) ; print(VP) u_850_1 = u850anom(:,{latS:latN},:) u_850 = u_850_1(lon|:,time|:, lat|:) U850 = dim_avg_Wrap(u_850(:,:,:)) u_200_1 = u200anom(lon|:,time|:, lat|:) u_200 = u_200_1(:,:,{latS:latN}) U200 = dim_avg_Wrap(u_200(:,:,:)) ; (time,lon) dimw = dimsizes( u200total ) ntim = dimw(0) nlat = dimw(1) mlon = dimw(2) ;delete(u200) f = addfile ("/ct13/ventrice/real_time/gfs200.nc" , "r") lon = f->lon time = f->time date = ut_calendar(time, -5) ; yyyymmddhh ;************************************************ ; Compute the temporal variance ;************************************************ var_vp = dim_variance_Wrap( VP) ; (lon) var_u850 = dim_variance_Wrap(U850) var_u200 = dim_variance_Wrap(U200) ;************************************************ ; Compute the zonal mean of the temporal variance ;************************************************ zavg_var_vp = dim_avg_Wrap( var_vp ) zavg_var_u850 = dim_avg_Wrap( var_u850) zavg_var_u200 = dim_avg_Wrap( var_u200) ;************************************************ ; Normalize by sqrt(avg_var*) ;************************************************ VP = VP/sqrt(zavg_var_vp ) ; (lon,time) U850 = U850/sqrt(zavg_var_u850) U200 = U200/sqrt(zavg_var_u200) ;************************************************ ; Combine the normalized data into one variable ;************************************************ printVarSummary (U850) printVarSummary (U200) printVarSummary(VP) cdata = new ( (/3*mlon,ntim/), typeof(VP), getFillValue(VP)) printVarSummary(cdata) do ml=0,mlon-1 cdata(ml ,:) = (/ VP(ml,:) /) cdata(ml+ mlon,:) = (/ U850(ml,:) /) cdata(ml+2*mlon,:) = (/ U200(ml,:) /) end do ;************************************************ ; Compute Combined EOF ;************************************************ ; eof_cdata = eofunc_Wrap(cdata , neof, False) ; (neof,3*mlon) fileName = "/ct13/ventrice/real_time/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 = 3 ; "vp", "u850", "u200" ceof = new( (/nvar,neof,mlon/), typeof(cdata), getFillValue(cdata)) do n=0,neof-1 ceof(0,n,:) = eof_cdata(n,0:mlon-1) ; vp ceof(1,n,:) = eof_cdata(n,mlon:2*mlon-1) ; u850 ceof(2,n,:) = eof_cdata(n,2*mlon:) ; u200 end do ceof_ts = new( (/nvar,neof,ntim/), typeof(cdata), getFillValue(cdata)) ceof_ts(0,:,:) = eofunc_ts_Wrap( VP,ceof(0,:,:),False) ; (neof,time) ceof_ts(1,:,:) = eofunc_ts_Wrap(U850,ceof(1,:,:),False) ; (neof,time) ceof_ts(2,:,:) = eofunc_ts_Wrap(U200,ceof(2,:,:),False) ; (neof,time) ;************************************************ ; Compute cross correlation of each variable's EOF time series at zero-lag ;************************************************ r_vp_u850 = escorc(ceof_ts(0,:,:) , ceof_ts(1,:,:)) r_vp_u200 = escorc(ceof_ts(0,:,:) , ceof_ts(2,:,:) ) r_u850_u200 = escorc(ceof_ts(1,:,:) , ceof_ts(2,:,:) ) print("==============") do n=0,neof-1 print("neof="+n \ +" r_vp_u850=" +sprintf("%4.3f",r_vp_u850(n)) \ +" r_vp_u200=" +sprintf("%4.3f",r_vp_u200(n)) \ +" r_u850_u200="+sprintf("%4.3f",r_u850_u200(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(1,:), mxlag) ; (N,mxlag+1) rlag_10 = esccr(eof_ts_cdata(1,:),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 = "/ct13/ventrice/real_time/" filo = "KELVIN_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 (pltType.eq."png") then pltTypeLocal = "eps" else pltTypeLocal = pltType 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" rts@xyExplicitLegendLabels = (/"VP", "U850", "U200" /) 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