Hello
Attached, please find an upodated
diagnostics_cam.ncl
[1] If this library is on your machine, then you can
replace it with the attached code.
[2] If you are not able to accomplish [1], then
extract the mjo_wavenum_freq_season_plot procedure
from the attached diagnostics_cam.ncl. Rename it to (say)
mjo_wavenum_freq_season_plot_traute,
then
load "./mjo_wavenum_freq_season_plot_traute"
and use this.
Questions, let me know.
Regards
Dennis Shea
Traute Crueger wrote:
> Hallo,
>
> I have a question with respect to the NCL MJO analysis tool (MJO Clivar).
> I applied the mjoclivar_10.ncl (wavenumber-frequency spectrum).
> The plot is perfect. However, in order to compare several of these
> plots, it would be useful that the colors levels are determined manually and not automatically.
> I made a few attemps, e.g. with
>
> optPlot_at_cnLevelSelectionMode="ManualLevels"
> optPlot_at_cnMinLevelValF = 0.0
> optPlot_at_cnMaxLevelValF = 8.0
> optPlot_at_cnLevelSpacingF = 0.5
> mjo_wavenum_freq_season_plot
> (wf,nameSeason(ns),pltDir,pltType,pltName,optPlot)
>
>
> however, the plot remains unchanged. How can I to solve this problem?
>
> Thank you and best regards
>
> Traute Crueger
>
;-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
; diagnostics_cam.ncl
;-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
; gsn_code.ncl, gsn_csm.ncl, contributed.ncl must be preloaded.
;-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
;
;-------------------------------------------------------------
; Decompose a variable into symmetric and asymmetric parts
;-------------------------------------------------------------
undef("decompose2SymAsym")
function decompose2SymAsym( z, iret[1] )
; supports: z(lat,lon), z(time,lat,lon), z(time,lev,lat,lon)
; antisymmetric part is stored in one hemisphere [eg: Northern Hemisphere]
; xOut( lat) = (x(lat)-x(-lat))/2
; symmetric part is stored in other hemisphere [eg: Southern Hemisphere]
; xOut(-lat) = (x(lat)+x(-lat))/2
local dimz, rankz, nlat, nl, N2, N, nt, kl, ntim, klev, dec2SymAsym, PRTFLG
begin
undef("dec2SymAsym") ; private to this bloc
function dec2SymAsym( x[*][*], nlat[1], N2[1], N[1] ) ; x(lat,lon)
; ===> local to decompose2SymAsym [~ f90 'contains' w PRIVATE attribute]
; This does the decomposition: called from driver 'decompose2SymAsym'.
local xsa, nl
begin
xsa = (/ x /) ; x(lat,lon)
do nl=0,N2-1
xsa( nl,:) = 0.5*(x(nlat-1-nl,:) + x(nl,:)) ; SH => symmetric
xsa(nlat-1-nl,:) = 0.5*(x(nlat-1-nl,:) - x(nl,:)) ; NH => asymmetric
end do
return( xsa )
end
; ====
dimz = dimsizes( z )
rankz = dimsizes( dimz )
if (rankz.ge.5) then
print("decompose2SymAsym: currently supports up to 4D: rank="+rankz+"D")
exit
end if
if (rankz.eq.1) then
nlat = dimz(0)
else
nlat = dimz(rankz-2)
end if
N2 = nlat/2
N = N2
if ((nlat%2).eq.1) then
N = N2+1 ; offset to handle the Equator
end if
zSymAsym = (/ z /) ; return array has same size/shape/type
; values will be overwritten
if (rankz.eq.1) then
do nl=0,N2-1
zSymAsym( nl) = 0.5*(z(nlat-1-nl) + z(nl)) ; SH => symmetric
zSymAsym(nlat-1-nl) = 0.5*(x(nlat-1-nl) - x(nl)) ; NH => asymmetric
end do
return( zSymAsym )
end if
if (rankz.eq.2) then
zSymAsym = dec2SymAsym( z, nlat, N2, N )
end if
if (rankz.eq.3) then
ntim = dimz(0)
do nt=0,ntim-1
zSymAsym(nt,:,:) = dec2SymAsym( z(nt,:,:), nlat, N2, N )
end do
end if
if (rankz.eq.4) then
ntim = dimz(0)
klev = dimz(1)
do nt=0,ntim-1
do kl=0,klev-1
zSymAsym(nt,kl,:,:) = dec2SymAsym( z(nt,kl,:,:), nlat, N2, N )
end do
end do
end if
zSymAsym_at_long_name = "antisymmetric & symmetric (separate hemispheres)"
if (isatt(z,"units")) then
zSymAsym_at_units = z_at_units
end if
copy_VarCoords( z, zSymAsym)
if (iret.eq.0)
return( zSymAsym ) ; original order
else
dNam = getvardims( zSymAsym )
if (rankz.eq.3) then
return( zSymAsym($dNam(1)$|:,$dNam(2)$|:,$dNam(0)$|:)) ; (lat,lon,time)
end if
if (rankz.eq.4) then
; return (lev,lat,lon,time)
return( zSymAsym($dNam(1)$|:,$dNam(2)$|:,$dNam(3)$|:,$dNam(0)$|:))
end if
end if
end
;-------------------------------------------------------------
; Prewhiten the data: eg remove the annual cycle.
; Actually, this will remove all time periods less than
; ththise corresponding to 'fCrit'.
; Note: The original fortran code provided by JET did not remove
; the grid point means so .... rmvMeans=False
; I assume that Matt Wheeler's code did that also.
;-------------------------------------------------------------
undef("rmvAnnualCycle")
function rmvAnnualCycle( z[*][*][*], spd[1], nDayTot[1]:integer
, rmvMeans[1]:logical \
, fCrit[1], iret[1] ) ; z(time,lat,lon)
local dimz, ntim, dNam, zt, cf, xbar
begin
dimz = dimsizes( z )
ntim = dimz(0)
dNam = getvardims( z ) ; input dimension names
zt = z($dNam(1)$|:,$dNam(2)$|:,$dNam(0)$|:) ; reorder: time fastest varying
cf = ezfftf( zt )
cf!0 = "cmplx" ; dimension clarity
cf!1 = "lat"
cf!2 = "lon"
cf!3 = "freq"
cf&cmplx = (/ 0,1 /)
cf&lat = z&$dNam(1)$
cf&lon = z&$dNam(2)$
; create freq dim values
freq = fspan(1,nDayTot*spd/2,nDayTot*spd/2)/nDayTot
freq!0 = "freq" ; clarity
freq_at_units = "cycles/day"
cf&freq = freq ; assign values, needed for {..}
cf(:,:,:,{:fCrit}) = 0.0 ; cf at all freq < fcrit = 0.0
; eg < 3/365 [6/730]
if (rmvMeans) then ; remove means for each lat,lon
; cf_at_xbar = 0.0
xbar = cf_at_xbar ; explicit retrieve
xbar = 0.0 ; set all means = 0.0
cf_at_xbar = xbar ; reassign new values
end if
; reconstruct
zt = (/ ezfftb(cf,cf_at_xbar) /) ; (/ ... /) avoid meta warning msg
if (iret.eq.0) then ; return in original order
return( zt($dNam(0)$|:,$dNam(1)$|:,$dNam(2)$|:) ) ; (time,lat,lon)
else ; return in permuted order
return( zt ) ; (lat,lon,time)
end if
end
;-------------------------------------------------------------
; Special reordering to resolve the Progressive and Retrogressive waves
; Reference: Hayashi, Y.
; A Generalized Method of Resolving Disturbances into
; Progressive and Retrogressive Waves by Space and
; Fourier and TimeCross Spectral Analysis
; J. Meteor. Soc. Japan, 1971, 49: 125-128.
;-------------------------------------------------------------
undef("resolveWavesHayashi")
function resolveWavesHayashi ( varfft[*][*][*], nDayWin[1], spd[1] ) ; (2,mlon,nSampWin)
;
; Create array PEE(NL+1,NT+1) which contains the (real) power spectrum.
; all the following assume indexing starting with 0
; In this array, the negative wavenumbers will be from pn=0 to NL/2-1;
; The positive wavenumbers will be for pn=NL/2+1 to NL.
; Negative frequencies will be from pt=0 to NT/2-1
; Positive frequencies will be from pt=NT/2+1 to NT .
; Information about zonal mean will be for pn=NL/2 .
; Information about time mean will be for pt=NT/2 .
; Information about the Nyquist Frequency is at pt=0 and pt=NT
;
; In PEE, define the
; WESTWARD waves to be either +ve frequency
; and -ve wavenumber or -ve freq and +ve wavenumber.
; EASTWARD waves are either +ve freq and +ve wavenumber
; OR -ve freq and -ve wavenumber.
; Note that frequencies are returned from fftpack are ordered like so
; input_time_pos [ 0 1 2 3 4 5 6 7 ]
; ouput_fft_coef [mean 1/7 2/7 3/7 nyquist -3/7 -2/7 -1/7]
; mean,pos freq to nyq,neg freq hi to lo
;
; Rearrange the coef array to give you power array of freq and wave number east/west
; Note east/west wave number *NOT* eq to fft wavenumber see Hayashi '71
; Hence, NCL's 'cfftf_frq_reorder' can *not* be used.
;
; For ffts that return the coefficients as described above, here is the algorithm
; coeff array varfft(2,n,t) dimensioned (2,0:numlon-1,0:numtim-1)
; new space/time pee(2,pn,pt) dimensioned (2,0:numlon ,0:numtim )
;
; NOTE: one larger in both freq/space dims
; the initial index of 2 is for the real (indx 0) and imag (indx 1) parts of the array
;
;
; if | 0 <= pn <= numlon/2-1 then | numlon/2 <= n <= 1
; | 0 <= pt < numtim/2-1 | numtim/2 <= t <= numtim-1
;
; if | 0 <= pn <= numlon/2-1 then | numlon/2 <= n <= 1
; | numtime/2 <= pt <= numtim | 0 <= t <= numtim/2
;
; if | numlon/2 <= pn <= numlon then | 0 <= n <= numlon/2
; | 0 <= pt <= numtim/2 | numtim/2 <= t <= 0
;
; if | numlon/2 <= pn <= numlon then | 0 <= n <= numlon/2
; | numtim/2+1 <= pt <= numtim | numtim-1 <= t <= numtim/2
local dimvf, numlon, N, varspacetime, pee, wave, freq
begin
dimvf = dimsizes( varfft )
mlon = dimvf(1)
N = dimvf(2)
varspacetime = new((/2,mlon+1,N+1/),typeof(varfft),1e20)
varspacetime(:,:mlon/2-1,:N/2-1) = varfft(:,mlon/2:1,N/2:N-1)
varspacetime(:,:mlon/2-1,N/2:) = varfft(:,mlon/2:1,:N/2)
varspacetime(:, mlon/2: ,:N/2) = varfft(:,:mlon/2,N/2:0)
varspacetime(:, mlon/2: ,N/2+1:) = varfft(:,:mlon/2,N-1:N/2)
; Create the real power spectrum pee = sqrt(real^2+imag^2)^2
pee = new((/mlon+1,N+1/),typeof(varfft),1e20)
pee = (sqrt(varspacetime(0,:,:)^2+varspacetime(1,:,:)^2))^2
; add meta data for use upon return
pee!0 = "wave"
pee!1 = "freq"
; JET fortran code 'wavep1'
wave = ispan(-mlon/2,mlon/2,1)
wave!0 = "wave"
wave&wave = wave
; JET fortran code 'frqfftwin'
freq = fspan(-1*nDayWin*spd/2,nDayWin*spd/2,nDayWin*spd+1)/nDayWin
freq!0 = "freq"
freq&freq = freq
pee&wave = wave
pee&freq = freq
return( pee )
end
;---------------------------------------------------------------
; dispersion curves
;--------------------------------------------------------------
undef("genDispersionCurves")
procedure genDispersionCurves( nWaveType[1]:integer \
, nEquivDepth[1]:integer \
, nPlanetaryWave[1]:integer \
, rlat[1]:numeric \
, Ahe[*]:numeric \
, Afreq[*][*][*]:numeric \
, Apzwn[*][*][*]:numeric )
local pi, re, g, omega, U, Un, ll, Beta, maxwn, ww, ed, wn, he \
, T, L, s, k, kn, del, deif, n, i, eif, P, dps, R, Rdeg, fillval
begin
pi = 4.0*atan(1.0)
re = 6.37122e06 ; [m] average radius of earth
g = 9.80665 ; [m/s] gravity at 45 deg lat used by the WMO
omega = 7.292e-05 ; [1/s] earth's angular vel
U = 0.0
Un = 0.0 ; since Un = U*T/L
ll = 2.*pi*re*cos(abs(rlat))
Beta = 2.*omega*cos(abs(rlat))/re
maxwn = nPlanetaryWave
fillval = 1e20
do ww = 1, nWaveType ; wave type
do ed = 1, nEquivDepth ; equivalent depth
he = Ahe(ed-1)
T = 1./sqrt(Beta)*(g*he)^(0.25)
L = (g*he)^(0.25)/sqrt(Beta)
do wn = 1, nPlanetaryWave ; planetary wave number
s = -20.*(wn-1)*2./(nPlanetaryWave-1) + 20.
k = 2.*pi*s/ll
kn = k*L
; Anti-symmetric curves
if (ww.eq.1) then ; MRG wave
if (k.lt.0.) then
del = sqrt(1.+(4.*Beta)/(k^2*sqrt(g*he)))
deif = k*sqrt(g*he)*(0.5-0.5*del)
end if
if (k.eq.0.) then
deif = sqrt(sqrt(g*he)*Beta)
end if
if (k.gt.0.) then
deif = fillval
end if
end if
if (ww.eq.2) then ; n=0 IG wave
if (k.lt.0.) then
deif = fillval
end if
if (k.eq.0.) then
deif = sqrt(sqrt(g*he)*Beta)
end if
if (k.gt.0.) then
del = sqrt(1.+(4.0*Beta)/(k^2*sqrt(g*he)))
deif = k*sqrt(g*he)*(0.5+0.5*del)
end if
end if
if (ww.eq.3) then ; n=2 IG wave
n=2.
del = (Beta*sqrt(g*he))
deif = sqrt((2.*n+1.)*del + (g*he)*k^2)
; do some corrections to the above calculated frequency.......
do i = 1,5
deif = sqrt((2.*n+1.)*del + (g*he)*k^2 + g*he*Beta*k/deif)
end do
end if
; symmetric curves
if (ww.eq.4) then ; n=1 ER wave
n=1.
if (k.lt.0.) then
del = (Beta/sqrt(g*he))*(2.*n+1.)
deif = -Beta*k/(k^2 + del)
else
deif = fillval
end if
end if
if (ww.eq.5) then ; Kelvin wave
deif = k*sqrt(g*he)
end if
if (ww.eq.6) then ; n=1 IG wave
n=1.
del = (Beta*sqrt(g*he))
deif = sqrt((2.*n+1.)*del + (g*he)*k^2)
; do some corrections to the above calculated frequency.......
do i=1,5
deif = sqrt((2.*n+1.)*del + (g*he)*k^2 + g*he*Beta*k/deif)
end do
end if
eif = deif ; + k*U since U=0.0
P = 2.*pi/(eif*24.*60.*60.)
dps = deif/k
R = L
Rdeg = (180.*R)/(pi*6.37e6)
Apzwn(ww-1,ed-1,wn-1) = s
if (deif.ne.fillval) then
P = 2.*pi/(eif*24.*60.*60.)
Afreq(ww-1,ed-1,wn-1) = 1./P
else
Afreq(ww-1,ed-1,wn-1) = fillval
end if
end do
end do
end do
end
;-------------------------------------------------------------
; Graphics Utility: Add temporal info the freq[y]- wave[x] contour
; freq is cpd [cycles per day]
;-------------------------------------------------------------
undef("addHorVertLines")
procedure addHorVertLines(wks[1]:graphic, plot[1]:graphic \
,x1[*],x2[*],y1[*],y2[*],y3[*],y4[*] )
; freq [y] axis: Add horizontal lines that explicitly
; print time in days. This assumes the units
; of the freq axis are "cpd" [cycles per day]
local gsres, txres
begin
gsres = True
gsres_at_gsLineDashPattern = 1
gsn_polyline(wks, plot, x1,y1,gsres)
gsn_polyline(wks, plot, x1,y2,gsres)
gsn_polyline(wks, plot, x1,y3,gsres)
gsn_polyline(wks, plot, x2,y4,gsres)
txres = True
txres_at_txJust = "CenterLeft"
txres_at_txFontHeightF = 0.013
gsn_text(wks, plot, "3 days" ,-14.7,.345,txres)
gsn_text(wks, plot, "6 days" ,-14.7,.175,txres)
gsn_text(wks, plot, "30 days",-14.7,.045,txres)
end
;-----------------------------------------------------------
; Utility: Unix does not allow spaces in file names
; When 5.0.1 becomes available, delete this and
; use replaceSingleChar in contributed.ncl
;-----------------------------------------------------------
undef("replaceChars")
procedure replaceChars (s[1]:string, oldStr[1]:string, newStr[1]:string )
local cOld, cNew, c, i
begin
cOld = stringtochar( oldStr )
cNew = stringtochar( newStr )
c = stringtochar( s )
i = ind( c.eq.cOld(0) )
if (.not.(any(ismissing(i)))) then
c(i) = cNew(0)
s = chartostring( c )
end if
end
;-----------------------------------------------------------
; Graphics Utility: This info may be of use if users want
; to set there own plot limits
; Activated if opt=True and (opt_at_debug=True or opt_at_Fig_3_statInfo=True)
;-----------------------------------------------------------
undef("statAsymSym")
procedure statAsymSym (x[*][*]:numeric, name[1]:string )
local statx, work, nwork
begin
statx = new ( 7, typeof(x), 1e20)
statx(0) = avg(x) ; mean
statx(1) = stddev(x) ; std. deviation
work = ndtooned(x)
nwork = dimsizes(work)
qsort( work ) ; sort into ascending order
statx(2) = min ( work )
statx(3) = work( nwork/4 ) ; lower quartile
statx(4) = work( nwork/2 ) ; approx median
statx(5) = work( nwork*3/4) ; upper quartile
statx(6) = max ( work )
print(" ")
print(" ===> Statistics: "+name+" <===")
print(" Mean="+statx(0) )
print(" StdDev="+statx(1) )
print(" Min="+statx(2) )
print(" lowQuartile="+statx(3) )
print(" Median="+statx(4) )
print("highQuartile="+statx(5) )
print(" Max="+statx(6) )
print(" ")
end
;==================================================================
; ====> Create Wheeler-Kiladis Space-Time plots. <====
;==================================================================
;==================================================================
; Note_1: The full logitudinal domain is used.
; This means that every planetary
; wavenumber will be represented.
; Note_2: Tapering in time is done to make the variable periodic.
;
; The calculations are also only made for the latitudes
; between '-latBound' and 'latBound'.
;
;******************** REFERENCES *******************************
; Wheeler, M., G.N. Kiladis
; Convectively Coupled Equatorial Waves:
; Analysis of Clouds and Temperature in the
; Wavenumber-Frequency Domain
; J. Atmos. Sci., 1999, 56: 374-399.
;---
; Hayashi, Y.
; A Generalized Method of Resolving Disturbances into
; Progressive and Retrogressive Waves by Space and
; Fourier and TimeCross Spectral Analysis
; J. Meteor. Soc. Japan, 1971, 49: 125-128.
;==================================================================
undef ("wkSpaceTime")
procedure wkSpaceTime (x[*][*][*]:numeric \
,diro[1]:string \
,caseName[1]:string \
,varName[1]:string \
,latBound[1]:numeric \
,spd[1]:integer \
,nDayWin[1]:integer \
,nDaySkip[1]:integer \
,opt[1]:logical)
local latN, latS, lonL, lonR, spd, fCrit, tim_taper \
, lon_taper, pltType, debug, minwav4smth, maxfrq4plt \
, minfrq4plt, maxfrq4plt, minwav4plt, maxwav4plt \
, fillVal, nMsg, dimx, ntim, nlat, mlon, nDayTot \
, nSampTot, nSampWin, nSampSkip, nWindow, N, dNam, work\
, rmvMeans, xAS, q, peeAS, nl, ntStrt, ntLast, nw, nt \
, ml, psumanti, psumsym, wv, wkdir, caseName, pltFilTit\
, frqfftwin, wavep1, minfrq, maxfrq, tmFontHgtF, pltTit\
, tiFontHgtF, lbFontHgtF, txFontHgtF, res, freq \
, wavenumber, NWVN, dcres, txres, rlat, Ahe \
, nWaveType, nPlanetaryWave, nEquivDepth, Apzwn, Afreq \
, asym, sym, x1, x2, y1, y2, y3, y4, wks, plot \
, psumb, psumsym_nolog, psumanti_nolog, tt, smthlen \
, i, pt8cpd, spec, nCn, nExtra
begin
debug = False ; default
if (opt .and. isatt(opt, "debug") ) then
debug = opt_at_debug
end if
if (isatt(x,"_FillValue")) then ; Check for _FillValue .... not allowed
nMsg = num(ismissing(x))
if (nMsg.gt.0) then
print("nMsg="+nMsg+" User must preprocess to remove _FillValue")
print(" FFTs do not allow missing values!! ")
exit
end if
delete(x@_FillValue) ; avoid warning messages from fft
end if
if (debug) then
printVarSummary( x )
printMinMax( x, True )
end if
;-------------------------------------------------------------------
; x sizes and dimension names
;-------------------------------------------------------------------
dimx = dimsizes(x)
ntim = dimx(0) ; total number of temporal samples
nlat = dimx(1)
mlon = dimx(2)
dNam = getvardims( x )
;-------------------------------------------------------------------
; check to make sure that "x" has full days of data
;-------------------------------------------------------------------
if ((ntim%spd).ne.0) then
nExtra = ntim%spd
print("nExtra="+nExtra+" input array must have complete days only ")
exit
end if
;-------------------------------------------------------------------
; Make input arguments into "internal" variables
;-------------------------------------------------------------------
latN = latBound
latS =-latBound ; make symmetric about the equator
lonL = 0 ; -180
lonR = 360 ; 180
fCrit = 1./nDayWin ; remove all contributions 'longer'
tim_taper = 0.1 ; time taper [0.1 => 10%]
lon_taper = 0.0 ; longitude taper [0.0 for globe; only global supported]
;-------------------------------------------------------------------
; Check for not allowed actions
;-------------------------------------------------------------------
if (lon_taper.gt.0.0 .or. (lonR-lonL).ne.360.) then
print("Code does currently allow lon_taper>0 or (lonR-lonL)<360")
exit
end if
;-------------------------------------------------------------------
; OPTIONS
;-------------------------------------------------------------------
pltType = "ps" ; default
if (opt .and. isatt(opt, "pltType") ) then
if (any(opt_at_pltType.eq.(/"ps", "eps", "x11", "ncgm", "png"/))) then
pltType = opt_at_pltType
end if
end if
pltTit = caseName+"_"+varName
pltTitle = pltTit+" LOG[Power: "+latBound+"S-"+latBound+"N]"
if (opt .and. isatt(opt, "pltTitle") ) then
pltTitle = opt_at_pltTitle
end if
pltFilTit = pltTit
replaceChars( pltFilTit, " ", "_") ; spaces not allowed unix file names
pltColorMap = "amwg_blueyellowred"
if (opt .and. isatt(opt, "pltColorMap") ) then
pltColorMap = opt_at_pltColorMap
end if
;-------------------------------------------------------------------
; Create required temporal sampling variables
;-------------------------------------------------------------------
nDayTot = ntim/spd ; # of days (total) for input variable
nSampTot = nDayTot*spd ; # of samples (total)
nSampWin = nDayWin*spd ; # of samples per temporal window
nSampSkip = nDaySkip*spd ; # of samples to skip between window segments
; neg means overlap
nWindow = (nSampTot-nSampWin)/(nSampWin+nSampSkip) + 1
N = nSampWin ; convenience [historical]
if (nDayTot.lt.nDayWin) then
print("nDayTot="+nDayTot+" is less the nDayWin="+nDayWin)
print(" This is not allowed !! ")
exit
end if
;-------------------------------------------------------------------
; Remove dominant signals
; (a) Explicitly remove *long term* linear trend
; For consistency with JET code keep the grid point means.
; This necessitates that 'dtrend_msg' be used because 'dtrend'
; always removes the mean(s).
; (b) All variations >= approx 'nDayWin' days if full year available
;-------------------------------------------------------------------
dNam = getvardims( x )
work = x($dNam(1)$|:,$dNam(2)$|:,$dNam(0)$|:) ; reorder (lat,lon,time)
;;work = dtrend( work , False ) ; remove mean + overall long term temporal trend
work = dtrend_msg(ispan(0,ntim-1,1)*1.0, work, False, False) ; remove just trend
if (isatt(work,"_FillValue")) then
delete(work@_FillValue) ; dtrend_msg adds this att
end if
; replace with detrended
x = (/ work($dNam(0)$|:,$dNam(1)$|:,$dNam(2)$|:) /) ; values (time,lat,lon)
delete(work)
if (nDayTot.ge.365) then ; rmv dominant signals
rmvMeans = False ; original code did not remove
x = rmvAnnualCycle(x, spd, nDayTot, rmvMeans, fCrit, 0)
end if
if (debug) then
print("===> Post removal of trend and signal <===")
printVarSummary( x ) ; (time,lat,lon)
printMinMax( x, True )
end if
;-------------------------------------------------------------------
; Decompose to Symmetric and Asymmetric parts
;-------------------------------------------------------------------
xAS = decompose2SymAsym( x , 1 ) ; create Asym and Sym parts
if (debug) then
printVarSummary(xAS) ; xAS(lat,lon,time) [iret=1]
printMinMax(xAS, True)
end if
;-------------------------------------------------------------------
; Because there is the possibility of overlapping *temporal* segments,
; we must use a less efficient approach and detrend/taper
; each window segment as it arises.
; t0 t1 t2 t3 t4 .................. t(N)
; lon(0): x00 x01 x02 x03 x04 .................. x0(N)
; : : : : : : :
; lon(M): xM0 xM1 xM2 xM3 xM4 .................. xM(N)
;-------------------------------------------------------------------
; q - temporary array to hold the 2D complex results
; for each longitude/time (lon,time) window that is fft'd.
; This is one instance [realization] of space-time decomposition.
;
; peeAS - symmetric and asymmetric power values in each latitude hemisphere.
; Add extra lon/time to match JET
;-------------------------------------------------------------------
q = new((/2,mlon,nSampWin/) ,typeof(xAS), "No_FillValue")
peeAS = new((/nlat,mlon+1,nSampWin+1/),typeof(xAS), "No_FillValue")
peeAS = 0.0 ; initialize
;-------------------------------------------------------------------
; Operate on the xAS array
; NCL does not have a complex 2D FFT at this time.
; Perform a "poorman's" complex 2D FFT by looping over time and space.
; Loop over all latitudes and then perform summing/averaging
; on the spectral results.
;-------------------------------------------------------------------
do nl=0,nlat-1
ntStrt = 0
ntLast = nSampWin-1
if (debug) then
print("==============> nl="+nl+" <==============")
end if
do nw=0,nWindow-1 ; temporal window
if (debug .and. nl.eq.0) then ; debug
print("nw="+nw+" ntStrt="+ntStrt+" ntLast="+ntLast)
end if
work = dtrend( xAS(:,:,ntStrt:ntLast), False ) ; detrend temporal window
work = taper ( work, tim_taper, 0) ; taper in "time" [periodic]
; work(nlat,mlon,N)
do nt=0,nSampWin-1 ; for each time perform complex fft in longitude
q(:,:,nt) = cfftf( work(nl,:,nt), 0.0, 0) ; space
end do
q = q/mlon ; normalize by # lon samples
do ml=0,mlon-1 ; for each lon perform complex fft in time
q(:,ml,:) = cfftf( q(0,ml,:), q(1,ml,:), 0) ; time
end do
q = q/nSampWin ; normalize by # time samples
;-------------------------------------------------------------------
; At this point 'q(2,mlon,nSampWin)' contains the
; real and imaginary space-time spectrum for this latitude
; ---
; Use Hayashi method to resolve into Progressive [Eastward]
; and Retrogressive [Westward] Waves.
;-------------------------------------------------------------------
pee = resolveWavesHayashi( q, nDayWin, spd )
peeAS(nl,:,:) = peeAS(nl,:,:) + (pee/nWindow) ; sum window contribution
ntStrt = ntLast+nSampSkip+1 ; set index for next temporal window
ntLast = ntStrt+nSampWin-1
end do ; nw (windows)
end do ; nl (lat)
delete(work)
;-------------------------------------------------------------------
; Add meta data to the Hayashi space-time symmetric and asymmetric power
;-------------------------------------------------------------------
peeAS!0 = "lat"
peeAS!1 = "wave"
peeAS!2 = "freq"
peeAS&lat = xAS&$dNam(1)$ ; nominally, xAS&lat
peeAS&wave = pee&wave
peeAS&freq = pee&freq
peeAS_at_long_name = "symmetric and asymmetric components"
if (debug) then
printVarSummary( peeAS )
printMinMax( peeAS , True)
end if
delete( q ) ; no longer needed
delete( pee )
;-------------------------------------------------------------------
; now that we have the power array for sym and asym: use to
; 1) plot raw power spectrum (some smoothing)
; 2) derive and plot the background spectrum (lots of smoothing)
; 3) derive a denoised spectrum that is raw power/background power
;-------------------------------------------------------------------
; psumanti and psumsym will contain the symmetric and asymmetric power
; summed over latitude
;-------------------------------------------------------------------
psumanti = new((/dimsizes(peeAS&wave),dimsizes(peeAS&freq)/),typeof(peeAS), 1e20 )
psumanti!0 = "wave"
psumanti!1 = "freq"
psumanti&wave = peeAS&wave
psumanti&freq = peeAS&freq
psumsym = psumanti
psumanti_at_long_name = "Asymmetric summed over lat"
psumsym_at_long_name = "Symmetric summed over lat"
if (nlat%2.eq.0) ; use named dimensions to reorder
psumanti = dim_sum_Wrap(peeAS(wave|:,freq|:,lat|nlat/2:nlat-1))
psumsym = dim_sum_Wrap(peeAS(wave|:,freq|:,lat|0:nlat/2-1) )
else
psumanti = dim_sum_Wrap(peeAS(wave|:,freq|:,lat|nlat/2+1:nlat-1))
psumsym = dim_sum_Wrap(peeAS(wave|:,freq|:,lat|0:nlat/2) )
end if
;-------------------------------------------------------------------
; since summing over half the array (symmetric,asymmetric) the
; total variance is 2x the half sum
;-------------------------------------------------------------------
psumanti = 2.0*psumanti
psumsym = 2.0*psumsym
;-------------------------------------------------------------------
; set the mean to missing to match original code
;-------------------------------------------------------------------
;
psumanti(:,{0.0}) = (/ psumanti@_FillValue /)
psumsym (:,{0.0}) = (/ psumsym@_FillValue /)
if (debug) then
printVarSummary( psumanti ) ; (wave,freq)
printMinMax( psumsym , True)
end if
;-------------------------------------------------------------------
; Apply smoothing to the spectrum. smooth over limited wave numbers
; Smoothing in frequency only (check if mean should be smoothed not smoothing now)
;--
; Smoothing parameters set these larger than the plotting
; wavenumbers to avoid smoothing artifacts
;-------------------------------------------------------------------
minwav4smth = -27
maxwav4smth = 27
do wv=minwav4smth,maxwav4smth
wk_smooth121( psumanti({wv},N/2+1:N-1) )
wk_smooth121( psumsym ({wv},N/2+1:N-1) )
end do
;-------------------------------------------------------------------
; Log10 scaling
;-------------------------------------------------------------------
psumanti_nolog = psumanti
psumsym_nolog = psumsym
psumanti = log10(psumanti)
psumsym = log10(psumsym)
psumanti_at_long_name = "log10(Asymmetric)"
psumsym_at_long_name = "log10(Symmetric)"
if (debug) then
printVarSummary( psumanti ) ; (wave,freq)
printMinMax( psumanti , True)
printVarSummary( psumsym )
printMinMax( psumsym , True)
end if
;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
; PLOT CODE follows
; --- set some 'plot variables
; set frequency maximum for plotting
; min(user specified frq, maxfrq in window)
;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
; The following allow DJS variable naming
; to be used with original plot code.
;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
;caseName = case
wkdir = diro
frqfftwin = peeAS&freq
frqfftwin&freq = peeAS&freq
wavep1 = peeAS&wave
wavep1&wave = peeAS&wave
if (debug) then
printVarSummary( frqfftwin )
printMinMax( frqfftwin, True )
printVarSummary( wavep1 )
printMinMax( wavep1, True )
end if
;-------------------------------------------------------------------
; plotting parameters freq and wavenumbers to plot
;-------------------------------------------------------------------
minfrq4plt = 0.
maxfrq4plt = 0.8
minwav4plt = -15
maxwav4plt = 15
minfrq = minfrq4plt
maxfrq = min((/maxfrq4plt,max(frqfftwin)/))
fillVal = 1e20 ; miscellaneous
;=============================================================
; Start Common Graphics Resources
;=============================================================
tmFontHgtF = 0.015 ; not sure why
tiFontHgtF = 0.018
lbFontHgtF = 0.015
txFontHgtF = 0.013
res = True
res_at_gsnFrame = False
res_at_gsnMaximize = True
res_at_gsnPaperOrientation = "portrait"
res_at_gsnLeftString = "Westward"
res_at_gsnRightString = "Eastward"
;res_at_lbBoxMinorExtentF = 0.18
res_at_lbLabelFontHeightF= lbFontHgtF
res_at_lbOrientation = "vertical"
res_at_cnFillOn = True
if (opt .and. isatt(opt, "cnLinesOn") \
.and. .not.opt_at_cnLinesOn) then
res_at_cnLinesOn = False
else
res_at_cnLineThicknessF = 0.5
end if
res_at_tmYLMode = "Explicit"
res_at_tmYLValues = fspan(minfrq,maxfrq,9)
res_at_tmYLLabels = fspan(minfrq,maxfrq,9)
res_at_tmYLMinorValues = fspan(minfrq,maxfrq,17)
res_at_tmYLLabelFontHeightF = tmFontHgtF
res_at_tmXBLabelFontHeightF = tmFontHgtF
res_at_tiXAxisString = "Zonal Wave Number"
res_at_tiXAxisFontHeightF= tiFontHgtF
res_at_tiYAxisString = "Frequency (cpd)"
res_at_tiYAxisFontHeightF= res_at_tiXAxisFontHeightF
if (.not.(pltTitle.eq."" .or. pltTitle.eq." ")) then
res_at_tiMainString = pltTitle
res_at_tiMainFontHeightF = tiFontHgtF
end if
res_at_txFontHeightF = tiFontHgtF
;------------------------------------------------------
; Create a list of variable names that have predefined
; contour intervals.
;------------------------------------------------------
varCnLevels = (/"FLUT" ,"OLR", "olr","U200","U850" \
,"PRECT","OMEGA500" /)
if (any(varCnLevels.eq.varName)) then
res_at_cnLevelSelectionMode = "ExplicitLevels"
else
res_at_cnLevelSelectionMode = "AutomaticLevels"
end if
;-------------------------------
; horizontal dashed lines and text for frequency axis [plot only]
;-------------------------------
freq = frqfftwin({freq|minfrq:maxfrq})
wavenumber = wavep1({wave|minwav4plt:maxwav4plt})
NWVN = dimsizes(wavenumber) ; number of wavenumbers
x1 = wavenumber
x1!0 = "wave"
y1 = new (NWVN,float)
y1!0 = "freq"
y1 = 1./3. ; 3 days
y2 = y1
y2 = 1./6. ; 6 days
y3 = y1
y3 = 1./30. ; 30 days
x2 = new(25,float)
x2 = 0.0
y4 = fspan(0.0,1.0,25)
;---------------------------------------------------------------
; dispersion: curves
;---------------------------------------------------------------
rlat = 0.0
Ahe = (/50.,25.,12./)
nWaveType = 6
nPlanetaryWave = 50
nEquivDepth = dimsizes(Ahe)
Apzwn = new((/nWaveType,nEquivDepth,nPlanetaryWave/),"double",fillVal)
Afreq = Apzwn
genDispersionCurves(nWaveType, nEquivDepth, nPlanetaryWave, rlat, Ahe, Afreq, Apzwn )
;---------------------------------------------------------------
; dispersion curve and text plot resources
;---------------------------------------------------------------
dcres = True
dcres_at_gsLineThicknessF = 2.0
dcres_at_gsLineDashPattern = 0
txres = True
txres_at_txJust = "CenterLeft"
txres_at_txPerimOn = True
txres_at_txFontHeightF = 0.013
txres_at_txBackgroundFillColor = "Background"
;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
; plotting params for fig 1; subset for plot
;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
; reorder so freq is "y"
asym = psumanti({freq|minfrq:maxfrq},{wave|minwav4plt:maxwav4plt})
asym!0 = "freq"
asym&freq = freq
asym!1 = "wave"
asym&wave = wavenumber
asym_at_long_name = "Fig_1: log10(Asymmetric)"
sym = psumsym({freq|minfrq:maxfrq},{wave|minwav4plt:maxwav4plt})
sym_at_long_name = "Fig_1: log10(Symmetric)"
if (debug) then
printVarSummary(asym)
printMinMax(asym, True)
printVarSummary(sym)
printMinMax(sym, True)
end if
;------------------------------------------------------
; Fig 1: Pre-defined contour levels [15] for selected variables [non-log10]
;------------------------------------------------------
nCn = 15
if (varName .eq. "FLUT" .or. varName .eq. "OLR" .or. varName .eq. "olr") then
res_at_cnLevels = (/-1.2,-1.1,-1.0,-0.8,-0.6,-0.4,-0.2 \ ; unequal
, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.1,1.2/)
end if
if (varName .eq. "PRECT") then
res_at_cnLevels = (/-18.2,-18.0,-17.8,-17.6,-17.5,-17.4,-17.3 \ ; unequal
,-17.2,-17.1,-17.0,-16.9,-16.8,-16.7,-16.6,-16.5/)
end if
if (varName .eq. "U200") then
res_at_cnLevels = fspan(-3.3, 0.9, nCn)
end if
if (varName .eq. "U850") then
res_at_cnLevels = fspan(-3.25, 0.25, nCn)
end if
if (varName .eq. "OMEGA500") then
res_at_cnLevels = fspan(-5.9, -4.5, nCn)
end if
if (opt .and. isatt(opt, "Fig_1")) then
res_at_cnLevelSelectionMode = "ExplicitLevels"
res_at_cnLevels = opt_at_Fig_1 ; user specified limits
end if
;------------------------------------------------------
; Fig 1: ANTI-SYMMETRIC
;------------------------------------------------------
;print("======> Fig 1a: ASYMMETRIC <=====")
wks = gsn_open_wks(pltType,wkdir+"/"+"Fig1.Asym."+pltFilTit)
gsn_define_colormap(wks,pltColorMap)
res_at_gsnCenterString = "Anti-Symmetric"
plot = gsn_csm_contour(wks,asym,res)
addHorVertLines(wks, plot, x1,x2,y1,y2,y3,y4 ) ; add dashed lines
frame(wks)
delete(wks) ; not required
;------------------------------------------------------
; Fig 1: SYMMETRIC
;------------------------------------------------------
;print("======> Fig 1b: SYMMETRIC <=====")
wks = gsn_open_wks(pltType,wkdir+"/"+"Fig1.Sym."+pltFilTit)
gsn_define_colormap(wks,pltColorMap)
res_at_gsnCenterString = "Symmetric"
plot = gsn_csm_contour(wks,sym,res)
addHorVertLines(wks, plot, x1,x2,y1,y2,y3,y4 ) ; add dashed lines
frame(wks)
delete(wks)
;------------------------------------------------------
; Is netCDF option set?
;------------------------------------------------------
if (opt .and. isatt(opt,"netCDF") .and. opt_at_netCDF) then
if (isatt(opt,"dirNetCDF")) then
dirNetCDF = opt_at_dirNetCDF
else
dirNetCDF = "./"
end if
if (isatt(opt,"filNetCDF")) then
filNetCDF = opt_at_filNetCDF
else
filNetCDF = "SpaceTime."+varName+".nc"
end if
fNam = dirNetCDF+filNetCDF
system ("/bin/rm -f "+fNam)
ncdf = addfile(fNam, "c")
ncdf->FIG_1_SYM = sym
ncdf->FIG_1_ASYM = asym
end if
;-----------------------------------------------------------------------------
; ****** now derive and plot the background spectrum (red noise) ************
; [1] Sum power over all latitude
; [2] Put fill value in mean
; [3] Apply smoothing to the spectrum. This smoothing DOES include wavenumber zero.
;-----------------------------------------------------------------------------
;print("======> BACKGROUND <=====")
psumb = dim_sum_Wrap(peeAS(wave|:,freq|:,lat|:)) ; sum over all latitudes
psumb_at_long_name = "Background Spectrum"
psumb@_FillValue = fillVal
psumb(wave|:,freq|N/2) = fillVal
psumb@_FillValue = fillVal
if (debug) then
printVarSummary(psumb) ; (wave,freq)
printMinMax(psumb, True)
end if
do tt = N/2+1,N
smthlen = maxwav4smth-minwav4smth+1
if (frqfftwin(tt).lt.0.1) then
do i = 1,5
wk_smooth121( psumb(freq|tt,{wave|minwav4smth:maxwav4smth}) )
end do
end if
if (frqfftwin(tt).ge.0.1.and.frqfftwin(tt).lt.0.2) then
do i = 1,10
wk_smooth121( psumb(freq|tt,{wave|minwav4smth:maxwav4smth}) )
end do
end if
if (frqfftwin(tt).ge.0.2.and.frqfftwin(tt).lt.0.3) then
do i = 1,20
wk_smooth121( psumb(freq|tt,{wave|minwav4smth:maxwav4smth}) )
end do
end if
if (frqfftwin(tt).ge.0.3) then
do i = 1,40
wk_smooth121(psumb(freq|tt,{wave|minwav4smth:maxwav4smth}))
end do
end if
end do
do nw = minwav4smth,maxwav4smth ; smth frequency up to .8 cycles per day
pt8cpd = min((/closest_val(.8,frqfftwin),dimsizes(frqfftwin)-1/))
smthlen = pt8cpd-(N/2+1)+1
do i = 1,10
wk_smooth121( psumb({nw},N/2+1:pt8cpd) )
end do
end do
;----------------------------------------------------------------
; [1] Put fill value in mean (again)
; [2] SAVE the background spectrum for plotting Fig 3
; [3] LOGARITHMIC SCALING for plotting the background spectrum
;----------------------------------------------------------------
psumb(wave|:,freq|N/2) = fillVal
psumb_nolog = psumb
psumb = log10(psumb) ; LOG10
psumb_at_long_name = "log10( background spec )"
if (debug) then
print(" ===> Post multiple smoothing by wk_smooth121 <===")
printVarSummary(psumb_nolog)
printMinMax(psumb_nolog, True)
printVarSummary(psumb)
printMinMax(psumb, True)
end if
;----------------------------------------------------------------
; set up for plotting [ subset of frequencies: not necessary]
;----------------------------------------------------------------
freq = frqfftwin({freq|minfrq:maxfrq})
spec = psumb({freq|minfrq:maxfrq},{wave|minwav4plt:maxwav4plt})
if (debug) then
printVarSummary(spec) ; (freq,wave)
printMinMax(spec, True)
end if
;----------------------------------------------------------------
; Fig 2: Predefined explicit contour levels BACKGROUND spectrum [LOG10]
;----------------------------------------------------------------
if (isatt(res,"cnLevels")) then
delete(res_at_cnLevels) ; allow size to change
end if
if (varName .eq. "FLUT" .or. varName .eq. "OLR" .or. varName .eq. "olr") then
res_at_cnLevels = (/-1.2,-1.1,-1.0,-0.8,-0.6,-0.4,-0.2 \ ; unequal 15
, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.1,1.2/)
end if
if (varName .eq. "PRECT") then
res_at_cnLevels = (/-18.2,-18.0,-17.8,-17.6,-17.5,-17.4,-17.3 \ ; unequal
,-17.2,-17.1,-17.0,-16.9,-16.8,-16.7,-16.6,-16.5/)
end if
if (varName .eq. "U200") then
res_at_cnLevels = fspan(-3.3, 0.9, nCn)
end if
if (varName .eq. "U850") then
res_at_cnLevels = fspan(-3.25, 0.25, nCn)
end if
if (varName .eq. "OMEGA500") then
res_at_cnLevels = fspan(-5.9,-4.5, nCn)
end if
if (opt .and. isatt(opt, "Fig_2")) then
res_at_cnLevelSelectionMode = "ExplicitLevels"
res_at_cnLevels = opt_at_Fig_2 ; user specified
end if
wks = gsn_open_wks(pltType,wkdir+"/Fig2."+pltFilTit)
gsn_define_colormap(wks,pltColorMap)
res_at_gsnCenterString = "Background Power"
plot = gsn_csm_contour(wks,spec,res)
addHorVertLines(wks, plot, x1,x2,y1,y2,y3,y4 ) ; add dashed lines
frame(wks)
;*********************************************************************************
; Fig 3a, 3b: psum_nolog/psumb_nolog [ratio]
;********************************************************************************
psumsym_nolog = psumsym_nolog /psumb_nolog ; (wave,freq)
psumanti_nolog = psumanti_nolog/psumb_nolog
if (debug) then
printVarSummary(psumanti_nolog) ; (freq,wave)
printMinMax(psumanti_nolog, True)
printVarSummary(psumsym_nolog) ; (freq,wave)
printMinMax(psumsym_nolog, True)
end if
;-----------------------------------------------------------
; ANTI-SYMMETRIC: RATIO ( subset to plot )
;-----------------------------------------------------------
asym = psumanti_nolog({freq|minfrq:maxfrq},{wave|minwav4plt:maxwav4plt})
sym = psumsym_nolog({freq|minfrq:maxfrq},{wave|minwav4plt:maxwav4plt})
if (debug) then
printVarSummary(asym) ; (freq,wave)
printMinMax(asym, True)
printVarSummary(sym) ; (freq,wave)
printMinMax(sym, True)
end if
if (isatt(res,"cnLevels")) then
delete(res_at_cnLevels) ; allow size to change
end if
if (varName .eq. "FLUT" .or. varName .eq. "OLR".or. varName .eq. "olr") then
res_at_cnLevels = fspan(0.3, 1.7, nCn)
end if
if (varName .eq. "U200") then
res_at_cnLevels = fspan(0.4, 1.8, nCn)
end if
if (varName .eq. "U850") then
res_at_cnLevels = fspan(0.4, 1.8, nCn)
end if
if (varName .eq. "PRECT") then
res_at_cnLevels = (/0.6,0.7 ,0.8,0.9 ,1.0,1.1,1.15,1.2,1.25 \
,1.3,1.35,1.4,1.45,1.5,1.6/)
end if
if (varName .eq. "OMEGA500") then
res_at_cnLevels = (/0.6,0.7,0.8,0.9,1.0,1.1,1.15,1.2,1.25 \
,1.3,1.35,1.4,1.45,1.5,1.6/)
end if
if (opt .and. isatt(opt, "Fig_3a")) then
res_at_cnLevelSelectionMode = "ExplicitLevels"
res_at_cnLevels = opt_at_Fig_3a ; user specified
end if
wks = gsn_open_wks(pltType,wkdir+"/Fig3.Asym."+pltFilTit)
gsn_define_colormap(wks,pltColorMap)
res_at_gsnCenterString = "Anti-Symmetric/Background"
plot = gsn_csm_contour(wks,asym,res)
addHorVertLines(wks, plot, x1,x2,y1,y2,y3,y4 ) ; add dashed lines
; draw dispersion line
gsn_polyline(wks,plot,Apzwn(0,0,:),Afreq(0,0,:),dcres)
gsn_polyline(wks,plot,Apzwn(0,1,:),Afreq(0,1,:),dcres)
gsn_polyline(wks,plot,Apzwn(0,2,:),Afreq(0,2,:),dcres)
gsn_polyline(wks,plot,Apzwn(1,0,:),Afreq(1,0,:),dcres)
gsn_polyline(wks,plot,Apzwn(1,1,:),Afreq(1,1,:),dcres)
gsn_polyline(wks,plot,Apzwn(1,2,:),Afreq(1,2,:),dcres)
gsn_polyline(wks,plot,Apzwn(2,0,:),Afreq(2,0,:),dcres)
gsn_polyline(wks,plot,Apzwn(2,1,:),Afreq(2,1,:),dcres)
gsn_polyline(wks,plot,Apzwn(2,2,:),Afreq(2,2,:),dcres)
; dispersion labels
gsn_text(wks,plot,"MRG",-10.0,.15,txres)
gsn_text(wks,plot,"n=2 IG",-3.0,.58,txres)
gsn_text(wks,plot,"n=0 EIG",9.5,.50,txres)
gsn_text(wks,plot,"h=50",-10.0,.78,txres)
gsn_text(wks,plot,"h=25",-10.0,.63,txres)
gsn_text(wks,plot,"h=12",-10.0,.51,txres)
frame(wks)
delete(wks) ; not required
;------------------------------------------------------------------
; SYMMETRIC
;------------------------------------------------------------------
if (isatt(res,"cnLevels")) then
delete(res_at_cnLevels)
end if
if (varName .eq. "FLUT" .or. varName .eq. "OLR".or. varName .eq. "olr") then
res_at_cnLevels = (/.3,.4,.5,.6,.7,.8,.9,1.,1.1,1.2,1.4,1.7,2.,2.4,2.8/)
end if
if (varName .eq. "U200") then
res_at_cnLevels = (/.4,.6,.8,1.,1.2,1.3,1.4,1.5,1.6,1.7,1.8,2.,2.2,2.4,2.6/)
end if
if (varName .eq. "U850") then
res_at_cnLevels = (/.4,.6,.8,1.,1.2,1.3,1.4,1.5,1.6,1.7,1.8,2.,2.2,2.4,2.6/)
end if
if (varName .eq. "PRECT") then
res_at_cnLevels = (/.6,.7,.8,.9,1.,1.1,1.15,1.2,1.25,1.3,1.35,1.4,1.45,1.5,1.6/)
end if
if (varName .eq. "OMEGA500") then
res_at_cnLevels = (/.6,.7,.8,.9,1.,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2./)
end if
if (opt .and. isatt(opt, "Fig_3b")) then
res_at_cnLevelSelectionMode = "ExplicitLevels"
res_at_cnLevels = opt_at_Fig_3b ; user specified
end if
wks = gsn_open_wks(pltType, wkdir+"/Fig3.Sym."+pltFilTit)
gsn_define_colormap(wks ,pltColorMap)
res_at_gsnCenterString = "Symmetric/Background"
plot = gsn_csm_contour(wks,sym,res)
addHorVertLines(wks, plot, x1,x2,y1,y2,y3,y4 ) ; add dashed lines
gsn_polyline(wks,plot,Apzwn(3,0,:),Afreq(3,0,:),dcres)
gsn_polyline(wks,plot,Apzwn(3,1,:),Afreq(3,1,:),dcres)
gsn_polyline(wks,plot,Apzwn(3,2,:),Afreq(3,2,:),dcres)
gsn_polyline(wks,plot,Apzwn(4,0,:),Afreq(4,0,:),dcres)
gsn_polyline(wks,plot,Apzwn(4,1,:),Afreq(4,1,:),dcres)
gsn_polyline(wks,plot,Apzwn(4,2,:),Afreq(4,2,:),dcres)
gsn_polyline(wks,plot,Apzwn(5,0,:),Afreq(5,0,:),dcres)
gsn_polyline(wks,plot,Apzwn(5,1,:),Afreq(5,1,:),dcres)
gsn_polyline(wks,plot,Apzwn(5,2,:),Afreq(5,2,:),dcres)
gsn_text(wks,plot,"Kelvin",11.0,.40,txres)
gsn_text(wks,plot,"n=1 ER",-10.7,.07,txres)
gsn_text(wks,plot,"n=1 IG",-3.0,.45,txres)
gsn_text(wks,plot,"h=50",-14.0,.78,txres)
gsn_text(wks,plot,"h=25",-14.0,.60,txres)
gsn_text(wks,plot,"h=12",-14.0,.46,txres)
frame(wks)
if (debug .or. (opt .and. isatt(opt,"Fig_3_statInfo") \
.and. opt_at_Fig_3_statInfo) ) then
statAsymSym( asym, "Fig_3a" )
statAsymSym( sym, "Fig_3b" )
end if
;------------------------------------------------------
; Is netCDF option set?
;------------------------------------------------------
if (opt .and. isatt(opt,"netCDF") .and. opt_at_netCDF) then
ncdf->FIG_3_BACK = spec
ncdf->FIG_3_SYM = sym
ncdf->FIG_3_ASYM = asym
end if
if (debug) then
print("********** FINISHED: Space-Time *****************")
end if
end
;-------------------------------------------------------------
; CAM interface to wkSpaceTime
;-------------------------------------------------------------
undef("wkSpaceTime_cam")
procedure wkSpaceTime_cam \
( diri[1]:string \ ; input directory
, fili[*]:string \ ; input file
, diro[1]:string \ ; output directory [plots]
, caseName[1]:string \ ; case name
, varName[1]:string \ ; variable name
, latBound[1]:numeric \ ; NH lat bound [positive]
, spd[1]:numeric \ ; samples per day
, level[1]:numeric \ ; if 4D ... what level?
, nDayWin[1]:integer \ ; temporal window length [days]
, nDaySkip[1]:integer \ ; days between time windows
, opt[1]:logical \ ; future options
)
;
; ---- local declarations are *not* required ----
local latN, latS, lonL, lonR, SPD, tskip, debug, fillVal \
, nfil, f, rank, nMsg, x, dNam, flist, dNamx \
, TIME, pltId
begin
latN = latBound
latS =-latBound ; make symmetric about the equator
lonL = 0 ; -180
lonR = 360 ; 180
; OPTIONS
debug = False ; default
if (opt .and. isatt(opt, "debug") ) then
debug = opt_at_debug
end if
SPD = spd ; SPD *actual* "samples per day"
; *after* possible 'decimation'
tskip = 1 ; used to skip samples on input
if (opt .and. isatt(opt, "spdSkip") .and. opt_at_spdSkip.gt.1 ) then
SPD = spd/opt_at_spdSkip
tskip = max( (/1,SPD/) ) ; decimate
end if
fillVal = 1e20 ; miscellaneous
;-------------------------------------------------------------------
; Read data from file(s): maintain meta data and original dim names
;-------------------------------------------------------------------
nfil = dimsizes(fili)
f = addfile (diri+fili(0), "r")
rank = dimsizes( filevardimsizes(f,varName) )
if (rank.lt.3 .or. rank.gt.4) then
print("wkSpaceTime: only 3D and 4D supported: rank="+rank+"D")
exit
end if
if (nfil.eq.1) then ; SINGLE FILE
if (rank.eq.3) then
x = f->$varName$(::tskip,{latS:latN},{lonL:lonR})
else
x = f->$varName$(::tskip,{level},{latS:latN},{lonL:lonR})
end if
else ; MULTIPLE FILES
dNam = getfilevardims(f,varName)
; needed for *large*
setfileoption("nc","SuppressClose",False) ; number of files
flist = addfiles( diri+fili, "r")
; make TIME [tedious]
TIME = flist[:]->$dNam(0)$ ; values
if (isfilevaratt(flist[0], dNam(0) , "units") ) then
TIME_at_units = flist[0]->$dNam(0)$@units ; assign units attribute
end if
if (isfilevarcoord( flist[0], dNam(0), dNam(0) ) ) then
TIME!0 = dNam(0) ; name the dimension
TIME&$dNam(0)$ = TIME ; assign values [coord]
end if
if (rank.eq.3) then
x = flist[:]->$varName$(::tskip,{latS:latN},{lonL:lonR})
x!0 = dNam(0)
x!1 = dNam(1)
x!2 = dNam(2)
else
x = flist[:]->$varName$(::tskip,{level},{latS:latN},{lonL:lonR})
x!0 = dNam(0)
x!1 = dNam(2)
x!2 = dNam(3)
end if
; NOT required but create meta data
dNamx = getvardims(x)
x&$dNamx(0)$ = TIME ; assign coordinates
x&$dNamx(1)$ = flist[0]->$dNamx(1)$({latS:latN})
x&$dNamx(2)$ = flist[0]->$dNamx(2)$({lonL:lonR})
if (isfilevaratt( flist[0], varName, "long_name")) then
x_at_long_name = flist[0]->$varName$@long_name
end if
if (isfilevaratt( flist[0], varName, "units" )) then
x_at_units = flist[0]->$varName$@units
end if
end if
;-------------------------------------------------------------------
; Call the procedure that will do the calculations and plots.
;-------------------------------------------------------------------
wkSpaceTime (x, diro, caseName, varName \
,latBound, SPD, nDayWin, nDaySkip, opt )
end
;-------------------------------------------------------
; L2 error norms
;-------------------------------------------------------
undef("l2_norm_a")
function l2_norma( x[*][*][*][*]:numeric, wy[*], wz )
local dimx, dimwz, rankwz, ntim, l2a, nt, xzon, xzonc \
, dif2, wyc, wywz, wywzs
begin
dimx = dimsizes( x )
dimwz = dimsizes( wz )
rankwz = dimsizes( dimwz )
ntim = dimx(0)
l2a = new ( ntim, typeof(x), getFillValue(x) )
; constant all times
wyc = conform(x(0,:,:,:), wy, 1) ; (klev,nlat,mlon)
if (rankwz.eq.1) then ; constant weight
wywz = wyc*conform( x(0,:,:,:), wz, 0 ) ; (klev,nlat,mlon)
wywzs = sum(wywz) ; scalar
delete( wyc )
end if
do nt=0,ntim-1
xzon = dim_avg( x(nt,:,:,:) ) ; (klev,nlat)
xzonc = conform( x(nt,:,:,:), xzon, (/0,1/)) ; (klev,nlat,mlon)
dif2 = (x(nt,:,:,:)-xzonc)^2 ; (klev,nlat,mlon)
if (rankwz.eq.1) then
l2a(nt) = sqrt( sum(dif2*wywz)/wywzs ) ; (nt)
else ; variable wgt (nt)
wywz = wyc*wz(nt,:,:,:) ; (klev,nlat,mlon)
l2a(nt) = sqrt( sum(dif2*wywz)/sum(wywz) ) ; (nt)
end if
end do
l2a_at_long_name = "l2"
if (isatt(x, "units")) then
l2a_at_units = x_at_units
end if
return( l2a )
end
;-------------------------------------------------------
; L2 error norms
;-------------------------------------------------------
undef("l2_norm_b")
function l2_normb( x[*][*][*][*]:numeric, wy[*], wz )
local dimx, dimwz, rankwz, ntim, l2b, nt, xzon0, xzonc \
, dif2, wyc, wywz, wywz
begin
dimx = dimsizes( x )
dimwz = dimsizes( wz )
rankwz = dimsizes( dimwz )
; time=0
xzon0 = dim_avg( x(0,:,:,:) ) ; (klev,nlat)
xzon0c = conform( x(0,:,:,:), xzon0, (/0,1/)) ; (klev,nlat,mlon)
wyc = conform( x(0,:,:,:), wy, 1) ; (klev,nlat,mlon)
if (rankwz.eq.1) then
wywz = wyc*conform( x(0,:,:,:), wz, 0 ) ; (klev,nlat,mlon)
wywzs = sum(wywz) ; scalar
delete( wyc )
end if
ntim = dimx(0)
l2b = new ( ntim, typeof(x), getFillValue(x) )
do nt=0,ntim-1
dif2 = (x(nt,:,:,:)-xzon0c)^2 ; (klev,nlat,mlon)
if (rankwz.eq.1) then
l2a(nt) = sqrt( sum(dif2*wywz)/wywzs ) ; (nt)
else
wywz = wyc*wz(nt,:,:,:) ; (klev,nlat,mlon)
l2b(nt) = sqrt( sum(dif2*wywz)/sum(wywz) ) ; (nt)
end if
end do
l2b_at_long_name = "l2"
if (isatt(x, "units")) then
l2b_at_units = x_at_units
end if
l2b_at_info = "x(:,:,:,:)-dimavg(x(0,:,:,:))"
return( l2b )
end
;-------------------------------------------------------
; l2 norm written by Mark Taylor
; function || T || on 2D horizontal array
;-------------------------------------------------------
undef ("norml2")
function norml2 (varz[*][*]:numeric,gw[*]:numeric)
local s2, gs, varl
begin
s2 = dimsizes(varz)
gs = dimsizes(gw)
if ( s2(0) .ne. gs(0) ) then
print ("norml2: error: first dimension does not match weight dimension: " \
+ s2(0) + " " + gs(0))
end if
varl = ( gw # (varz^2) )/sum(gw)
return( sqrt( sum(varl)/s2(1) ))
end
;-------------------------------------------------------
; Energy
;-------------------------------------------------------
undef ("energy_a")
function energy_a(u[*][*][*][*]:numeric, v[*][*][*][*]:numeric\
,t[*][*][*][*]:numeric,phis[*][*][*]:numeric \
,wy[*]:numeric, wz:numeric, opt[1]:integer )
;-------------------------------------------------------
; energy = SUM(0.5*(U^2+V^2) + cp*T + S )*wz
;-------------------------------------------------------
local dimu, dimwz, rankwz, ntim, e, nt, x, wyc, wywz, wywzs
begin
dimu = dimsizes( u )
dimwz = dimsizes( wz )
rankwz = dimsizes( dimwz )
ntim = dimu(0)
e = new ( ntim, typeof(u), getFillValue(u) )
; constant all times
wyc = conform(u(0,:,:,:), wy, 1) ; (klev,nlat,mlon)
if (rankwz.eq.1) then ; constant weight
wywz = wyc*conform( u(0,:,:,:), wz, 0 ) ; (klev,nlat,mlon)
wywzs = sum(wywz) ; scalar
delete( wyc )
end if
a = 6.371229e6 ; m
g = 9.80616 ; m/s^2
cp = 1004.64 ; J/(kg-K)
twopi = 8.*atan(1.0)
do nt=0,ntim-1
s = conform(u(nt,:,:,:), phis(nt,:,:), (/1,2/))
x = (0.5*(u(nt,:,:,:)^2+v(nt,:,:,:)^2) + cp*t(nt,:,:,:) + s )
if (rankwz.eq.1) then
e(nt) = sum(x*wywz)
else ; variable wgt (nt)
wywz = wyc*wz(nt,:,:,:) ; (klev,nlat,mlon)
e(nt) = sum(x*wywz)
end if
end do
e = (twopi*a*a/g)*e ; total energy e(t)
if (opt.eq.0) then
e_at_long_name = "E: eqn 155"
e_at_units = "J" ; ==>'kg m^2 / s^2'.
else
e = (e-e(0))/e(0) ; normalized
e_at_long_name = "Normalized Energy Difference"
e_at_units = "(E-E(0))/E(0)"
end if
return( e )
end
; ====================================================
undef("mjo_ApplyLanczosWgt_LeftDim")
function mjo_ApplyLanczosWgt_LeftDim(x:numeric, wgt[*]:numeric, opt[1]:integer)
;
; utility routine ... makes for cleaner code
; Reorder so time is rightmost; applies Lanczos weights, reorder back
;
local dimx, rank, dNam, xBand
begin
dimx = dimsizes(x)
rank = dimsizes( dimx )
dNam = getvardims( x )
if (rank.eq.1) then
return( wgt_runave_Wrap( x, wgt, opt) )
end if
if (rank.eq.2) then
xBand = wgt_runave_Wrap( x($dNam(1)$|:,$dNam(0)$|:), wgt, opt)
return( xBand($dNam(0)$|:,$dNam(1)$|:) )
end if
if (rank.eq.3) then
xBand = wgt_runave_Wrap( x($dNam(1)$|:,$dNam(2)$|:,$dNam(0)$|:), wgt, opt)
return( xBand($dNam(0)$|:,$dNam(1)$|:,$dNam(2)$|:) )
end if
if (rank.eq.4) then
xBand = wgt_runave_Wrap( x($dNam(1)$|:,$dNam(2)$|:,$dNam(3)$|:,$dNam(0)$|:), wgt, opt)
return (xBand($dNam(0)$|:,$dNam(1)$|:,$dNam(2)$|:,$dNam(3)$|:) )
end if
if (rank.eq.5) then
xBand = wgt_runave_Wrap( x($dNam(1)$|:,$dNam(2)$|:,$dNam(3)$|:,$dNam(4)$|:,$dNam(0)$|:), wgt, opt)
return (xBand($dNam(0)$|:,$dNam(1)$|:,$dNam(2)$|:,$dNam(3)$|:,$dNam(4)$|:) )
end if
end
; ====================================================
undef("apply_dtrend_LeftDim")
function apply_dtrend_LeftDim(x:numeric, wgt[*]:numeric, opt[1]:integer)
;
; utility routine ... makes for cleaner code
; Reorder so time is rightmost; detrend the time series, reorder back
;
local dimx, rank, dNam, x_dtrend, xWork
begin
dimx = dimsizes(x)
rank = dimsizes( dimx )
dNam = getvardims( x )
if (rank.eq.1) then
x_dtrend = x ; retain meta data
x_dtrend = dtrend( x, False)
end if
if (rank.eq.2) then
xWork = x($dNam(1)$|:,$dNam(0)$|:) ; retain meta data
xWork = dtrend( xWork, False)
return( xWork($dNam(0)$|:,$dNam(1)$|:) )
end if
if (rank.eq.3) then
xWork = x($dNam(1)$|:,$dNam(2)$|:,$dNam(0)$|:)
xWork = dtrend( xWork, False)
return( xWork($dNam(0)$|:,$dNam(1)$|:,$dNam(2)$|:) )
end if
if (rank.eq.4) then
xWork = x($dNam(1)$|:,$dNam(2)$|:,$dNam(3)$|:,$dNam(0)$|:)
xWork = dtrend( xWork, False)
return (xWork($dNam(0)$|:,$dNam(1)$|:,$dNam(2)$|:,$dNam(3)$|:) )
end if
if (rank.eq.5) then
xWork = x($dNam(1)$|:,$dNam(2)$|:,$dNam(3)$|:,$dNam(4)$|:,$dNam(0)$|:)
xWork = dtrend( xWork, False)
return (xWork )
end if
end
; ====================================================
undef ("band_pass_area_time")
function band_pass_area_time \
(x[*][*][*]:numeric, spd[1]:numeric \
,bpf[3]:numeric \
,wy[*]:numeric, opt:logical)
local nMsg, dimx, dimwy, ntim, bpfStrt, bpfLast, bpfNwgt \
, dumy, save_FillValue, fca, fcb, ihp, sigma, xts
begin
; error check
nMsg = num( ismissing(x) )
if (nMsg.gt.0) then
print("band_pass_area_time: currently, missing data not allowed: nMsg="+nMsg)
exit
end if
dimx = dimsizes(x)
dimwy = dimsizes(wy)
if (dimx(1).ne.dimwy(0)) then
print("band_pass_area_time: sizes of y/lat dimension do not match")
print(" dimsizes(wy)="+dimwy)
print(" dimx(1)=nlat="+dimx(1))
exit
end if
dimx = dimsizes( x )
ntim = dimx(0)
bpfStrt = bpf(0) ; days
bpfLast = bpf(1)
bpfNwgt = bpf(2) ; effective # weights
bpfNwgt = (bpfNwgt/2)*2 + 1 ; make sure it is odd
if (bpfStrt.gt.bpfLast) then ; safety check
dumy = bpfLast
bpfLast = bpfStrt
bpfStrt = dumy
end if
fca = 1.0/(spd*bpfLast) ; freq start/end
fcb = 1.0/(spd*bpfStrt) ; bpf{Strt/Last}
if (typeof(x).eq."double" .or. typeof(wy).eq."double") then
xts = new( (/3,ntim/), "double", getFillValue(x))
else
xts = new( (/3,ntim/), "float" , getFillValue(x))
end if
xts!0 = "var"
; weighted area averages
xts(1,:) = wgt_areaave_Wrap(x, wy, 1., 0) ; (time)
if (opt .and. isatt(opt,"detrend") .and. opt_at_detrend) then
if (isatt(xts,"_FillValue")) then
save_FillValue = x@_FillValue
delete(xts@_FillValue) ; avoid annoying warning msg
end if ; from dtrend
xts(1,:) = dtrend(xts(1,:), False ) ; detrend
if (isvar("save_FillValue")) then
xts@_FillValue = save_FillValue ; reassign
end if
end if
ihp = 2 ; bpf=>band pass filter
sigma = 1.0 ; Lanczos sigma
bpfwgt = filwgts_lanczos (bpfNwgt, ihp, fca, fcb, sigma )
xts(0,:) = wgt_runave_Wrap(xts(1,:), bpfwgt, 0)
if (opt .and. isatt(opt,"Nrun")) then
Nrun = (opt_at_Nrun/2)*2 + 1 ; make sure it is odd
else
Nrun = 101
end if
len = (Nrun*spd)/2 ; half total length
do n = len,ntim-(len+1)
xts(2,n) = variance(xts(0,n-len:n+len)) ; 'local' variance
end do
xts_at_bpfStrt = bpfStrt
xts_at_bpfLast = bpfLast
xts_at_bpfNwgt = bpfNwgt
xts_at_var_0 = "band pass"
xts_at_var_1 = "raw areal means"
xts_at_var_2 = "local variances: Nrun="+Nrun
return( xts )
end
; =================================================================
undef("band_pass_area_time_plot")
procedure band_pass_area_time_plot( \
xts[3][*]:numeric, time[*]:numeric, \
pltDir[1]:string, pltType[1]:string,\
pltName[1]:string, opt[1]:logical)
begin
date = ut_calendar(time, -3)
yrfrac = yyyymmddhh2yyyyFrac( date)
delete(yrfrac_at_long_name)
pltPath= pltDir+"/"+pltName
wks = gsn_open_wks(pltType , pltPath)
gsn_merge_colormaps(wks,"amwg_blueyellowred","BlAqGrYeOrReVi200")
;gsn_draw_colormap(wks) ; draw color map
plot = new ( 3, "graphic")
; left variable
res = True
res_at_gsnDraw = False
res_at_gsnFrame = False
res_at_vpXF = 0.15
res_at_vpYF = 0.3
res_at_vpWidthF = 0.65
res_at_vpHeightF = 0.18
res_at_tmXBLabelsOn = False ; do not draw bottom labels
res_at_tmXBOn = False ; no bottom tickmarks
res_at_xyLineThicknessF = 1. ; 1 is default
res_at_tiYAxisString = "Raw areal Avg"
plot(0) = gsn_csm_xy(wks,yrfrac,xts(1,:),res)
res_at_gsnYRefLine = 0.
res_at_gsnAboveYRefLineColor = "Red"
res_at_gsnBelowYRefLineColor = "Blue"
res_at_tiYAxisString = "Band-Pass"
plot(1) = gsn_csm_xy(wks,yrfrac,xts(0,:),res)
delete(res_at_gsnAboveYRefLineColor)
delete(res_at_gsnBelowYRefLineColor)
delete(res_at_gsnYRefLine)
res_at_tmXBLabelsOn = True
res_at_tmXBOn = True
res_at_tmXBFormat = "f" ; remove trailing .0s
res_at_xyLineColor = "green"
res_at_xyLineThicknessF = 2. ; 1 is default
yrfrac_at_long_name = "date as fraction of year"
res_at_tiYAxisString = "Run Variance"
plot(2) = gsn_csm_xy(wks,yrfrac,xts(2,:),res)
;************************************************
resP = True ; modify the panel plot
resP_at_txString = xts_at_long_name+": Areal Averaged & Filtered: " \
+xts_at_bpfStrt+"-"+xts_at_bpfLast+" days: nw=" \
+xts_at_bpfNwgt
resP_at_gsnMaximize = True
resP_at_gsnPaperOrientation = "portrait" ; force portrait
resP_at_gsnPanelBottom = 0.05 ; add some space at bottom
gsn_panel(wks,plot,(/3,1/),resP)
end
; +++++
undef ("band_pass_area_time_plot_cam")
procedure band_pass_area_time_plot_cam( \
xts[3][*]:numeric, date[*]:integer, datesec[*]:integer, \
pltDir[1]:string, pltType[1]:string, pltName[1]:string)
begin
end
;=========================================
undef ("band_pass_latlon_time")
function band_pass_latlon_time \
(x[*][*][*]:numeric, spd[1]:numeric \
,bpf[3]:numeric, opt:logical)
local nMsg, dimx, ntim, nlat, mlon, dNam, n, tim_taper \
, bpfStrt, bpfLast, bpfNwgt, fca, fcb, ihp, sigma \
, bpfWgt, work, WORK, cf, WCF
begin
nMsg = num( ismissing(x) )
if (nMsg.gt.0) then
print("band_pass_latlon_time: currently, missing data not allowed: nMsg="+nMsg)
exit
end if
dimx = dimsizes(x)
ntim = dimx(0)
nlat = dimx(1)
mlon = dimx(2)
bpfStrt = bpf(0) ; days
bpfLast = bpf(1)
bpfNwgt = bpf(2) ; effective # weights
bpfNwgt = (bpfNwgt/2)*2 + 1 ; make sure it is odd
if (bpfStrt.gt.bpfLast) then ; safety check
dumy = bpfLast
bpfLast = bpfStrt
bpfStrt = dumy
end if
fca = 1.0/(spd*bpfLast)
fcb = 1.0/(spd*bpfStrt)
dNam = getvardims(x) ; get dimension names
do n=0,2 ; only used if detrending
if (ismissing(dNam(n))) then ; or tapering in time
x!n = "dim"+n ; assign name
dNam = x!n
end if
end do
;;work = x(lat|:,lon|:,time|:) ; reorder so time is rightmost
work = x($dNam(1)$|:, $dNam(2)$|:, $dNam(0)$|:) ; generic
; By default ... no detrend
if (opt .and.(isatt(opt,"detrend") .and. opt_at_detrend)) then
if (isatt(work,"_FillValue")) then
delete(work@_FillValue) ; avoid annoying warning msg
end if
work = dtrend(work, False ) ; detrend
work_at_detrend = "data detrended in time"
end if
ihp = 2 ; bpf=>band pass filter
sigma = 1.0 ; Lanczos sigma
bpfwgt = filwgts_lanczos (bpfNwgt, ihp, fca, fcb, sigma )
if (opt .and. isatt(opt,"fft") .and. opt_at_fft) then
; By default ... no taper
if (isatt(opt,"taper")) then
tim_taper = opt_at_taper
else
tim_taper = 0.10 ; default taper is 10%
end if
work = taper(work, tim_taper, 0)
work_at_taper = "data tapered in time"
; fft in time
cf = ezfftf (work) ; cf(2,nlat,mlon,ntim/2)
; map response of digitial
; filter to fft space
fcf = fspan(0, 0.5, ntim/2) ; fft freq
wcf = linint1 (bpfwgt_at_freq, bpfwgt_at_resp, False, fcf, 0)
WCF = conform(cf(0,:,:,:), wcf, 2)
delete(wcf)
cf(0,:,:,:) = cf(0,:,:,:)*WCF ; apply response coef
cf(1,:,:,:) = cf(1,:,:,:)*WCF
delete(WCF)
work = ezfftb(cf, 0.0) ; fourier synthesis
delete(cf)
work_at_process = "FFT with digital respoonse mapped tp FFT space"
else
work = wgt_runave_Wrap(work,bpfwgt,0)
work_at_process = "wgt_runave"
end if
work_at_band_pass_start = bpfStrt
work_at_band_pass_last = bpfLast
work_at_band_pass_Nwgts = bpfNwgt
X = work($dNam(0)$|:, $dNam(1)$|:, $dNam(2)$|:)
copy_VarMeta (x, X)
X_at_long_name = "Band Pass: "
if (isatt(x,"long_name")) then
X_at_long_name = "Band Pass: "+x_at_long_name
end if
return( X )
end
; ----------------------------------------------------
undef ("band_pass_hovmueller")
function band_pass_hovmueller( \
x[*][*][*]:numeric, spd[1]:numeric \
,bpf[3]:numeric, wy[*]:numeric, opt:logical)
local nMsg, dimx, ntim, nlat, mlon, dNam, n, saveFillV \
, bpfStrt, bpfLast, bpfNwgt, fca, fcb, ihp, sigma \
, bpfWgt, wgty, work, WORK, cf, WCF
begin
nMsg = num( ismissing(x) ) ; error check
if (nMsg.gt.0) then
print("band_pass_hovmueller: currently, missing data not allowed: nMsg="+nMsg)
exit
end if
dimx = dimsizes(x)
dimwy = dimsizes(wy) ; size wy
if (dimwy.gt.1 .and. dimx(1).ne.dimwy(0)) then
print("band_pass_hovmueller: sizes of y/lat dimension do not match")
print(" dimsizes(wy)="+dimwy)
print(" dimx(1) ="+dimx(1))
exit
end if
if (dimwy.eq.1) then ; scalar
wgty = new( dimwy, typeof(wy), getFillValue(wy))
end if
wgty = wy
ntim = dimx(0)
nlat = dimx(1)
mlon = dimx(2)
bpfStrt = bpf(0) ; days
bpfLast = bpf(1)
bpfNwgt = bpf(2) ; effective # weights
bpfNwgt = (bpfNwgt/2)*2 + 1 ; make sure it is odd
if (bpfStrt.gt.bpfLast) then ; safety check
dumy = bpfLast
bpfLast = bpfStrt
bpfStrt = dumy
end if
fca = 1.0/(spd*bpfLast) ; frq corresponding to
fcb = 1.0/(spd*bpfStrt) ; bpfStrt and bpfLast
dNam = getvardims(x) ; get dimension names
do n=0,2 ; only used if detrending
if (ismissing(dNam(n))) then
x!n = "dim"+n
dNam = x!n
end if
end do
; weighted average at each (lon,time)
work = dim_avg_wgt_Wrap(x($dNam(2)$|:, $dNam(0)$|:, $dNam(1)$|:), wgty, 0)
; unweighted average at each (lon,time)
;work = dim_avg_Wrap(x($dNam(2)$|:, $dNam(0)$|:, $dNam(1)$|:))
; Remove overall mean
work = work - avg(work)
; By default ... no detrend in time
if (opt .and.(isatt(opt,"detrend") .and. opt_at_detrend)) then
if (isatt(work,"_FillValue")) then
saveFillV = work@_FillValue ; avoid annoying warning msg
delete(work@_FillValue)
end if
work = dtrend(work, False ) ; detrend in time
if (isvar("saveFillV")) then
work@_FillValue = saveFillV ; reassign
end if
work_at_detrend = "data detrended in time"
end if
ihp = 2 ; bpf=>band pass filter
sigma = 1.0 ; Lanczos sigma
bpfwgt = filwgts_lanczos (bpfNwgt, ihp, fca, fcb, sigma )
work = wgt_runave_Wrap(work,bpfwgt,0) ; apply filter to time
work_at_band_pass_start = bpfStrt
work_at_band_pass_last = bpfLast
work_at_band_pass_Nwgts = bpfNwgt
return( work($dNam(0)$|:, $dNam(2)$|:) ) ; (time,lon)
end
; ==============> prototype plot procedure <=================
; create Hovmueller plots
; =============================================================
undef ("band_pass_hovmueller_plot")
procedure band_pass_hovmueller_plot( x[*][*]:numeric \ ; (time,lon)
, pltDir[1]:string \
, pltType[1]:string \
, pltName[1]:string \
, opt[1]:logical )
local date, yrfrac, TIME, ntim, yyyy, mm, dd, hh, pltPath \
, resH, hplt, nt, ntm, monNam, tValL, tLabL \
, tValR, tLabR
begin
pltPath = pltDir+"/"+pltName
wks = gsn_open_wks(pltType , pltPath)
gsn_define_colormap(wks,"amwg_blueyellowred")
resH = True ; plot mods desired
resH_at_cnFillOn = True ; turn on color fill
resH_at_cnLinesOn = False ; turn of contour lines
resH_at_cnLineLabelsOn = False ; turn of contour line labels
resH_at_gsnSpreadColors = True ; use full range of color map
resH_at_gsnSpreadColorStart = 2
resH_at_gsnSpreadColorEnd = 17
resH_at_lbLabelAutoStride = True ; let NCL figure spacing
resH_at_pmLabelBarHeightF = 0.1 ; default is taller
resH_at_trYReverse = True ; reverse the y (time) axis
resH_at_gsnMaximize = True
resH_at_gsnPaperOrientation = "portrait"
symMinMaxPlt (x,16,False,resH) ; nice symmetric range
resH_at_tmLabelAutoStride = True ; don't let labels overlap
; TIME RELATED VARIABLES
time = x&time ; clarity
date = ut_calendar(x&time, -3) ; gregorian calendar
yrfrac = yyyymmddhh2yyyyFrac( date)
delete(yrfrac_at_long_name)
TIME = ut_calendar(x&time, 0 ) ; gregorian calendar
yyyy = floattoint( TIME(:,0) )
mm = floattoint( TIME(:,1) )
dd = floattoint( TIME(:,2) )
;hh = floattoint( TIME(:,3) )
ntim = dimsizes(date)
; Left axis
if (opt .and. isatt(opt,"yearFraction")) then
resH_at_tmYLMode = "Manual"
if (isatt(opt,"yearFractionSpacingF")) then
resH_at_tmYLTickSpacingF = opt_at_yearFractionSpacingF
else
resH_at_tmYLTickSpacingF = 0.25 ; 0, .25, .50, .75
end if
;;resH_at_tmYLTickStartF = yyyy(0)
;;resH_at_tmYLTickEndF = yyyy(ntim-1)
x&time = yrfrac
else
monNam = (/"Jan","Feb","Mar","Apr","May","Jun" \
,"Jul","Aug","Sep","Oct","Nov","Dec" /)
tValL = new(ntim, typeof(x&time) , "No_FillValue") ; bigger than
tLabL = new(ntim, "string", "No_FillValue") ; needed
;tValR = tValL
;tLabR = tLabL
if (opt .and. isatt(opt,"monthLabels")) then
monthLabels = opt_at_monthLabels
else
monthLabels = (/1,4,7,10/)
end if
ntm = -1
do nt=0,ntim-1
if (dd(nt).eq.1 .and. any(mm(nt).eq.monthLabels)) then
ntm = ntm + 1
tValL(ntm) = x&time(nt)
tLabL(ntm) = monNam(mm(nt)-1)
;tValR(ntm) = x&time(nt)
if (dd(nt).eq.1 .and. mm(nt).eq.1) then
tLabL(ntm) = sprinti("%0.4i", yyyy(nt))+" "+tLabL(ntm)
; tLabL(ntm) = yyyy(nt)+" "+tLabL(ntm)
; ;tLabL(ntm) = tLabL(ntm)+"~C~"+yyyy(nt)
; tLabR(ntm) = yyyy(nt)
;else
; tLabR(ntm) = ""
end if
end if
end do
resH_at_tmYLMode = "Explicit"
resH_at_tmYLValues = tValL(0:ntm)
resH_at_tmYLLabels = tLabL(0:ntm)
;resH_at_tmYRLabelsOn = True ; use for year
;resH_at_tmYUseLeft = False
;resH_at_tmYRMode = "Explicit"
;resH_at_tmYRValues = tValR(0:ntm)
;resH_at_tmYRLabels = tLabR(0:ntm)
end if
;resH_at_cnLevelSpacingF = 20
;resH_at_cnMaxLevelValF = 80
;resH_at_cnMinLevelValF = -80
;resH_at_cnLevelSelectionMode = "ManualLevels"
;hplt = gsn_csm_hov(wks, x, resH)
hplt = gsn_csm_contour(wks, x, resH)
if (isatt(resH,"tmYLMode") .and. resH_at_tmYLMode.eq."Manual") then
x&time = time ; reassign original time coordinate variable
end if
end
; ====================================================
undef("mjo_spectra_season")
function mjo_spectra_season \
(x[*][*][*]:numeric, date[*]:integer, wy[*]:numeric \
,seaName[1]:string, opt:logical)
; MJO CLIVAR: perform segment averaging
; winter starts Nov 1 [180 days]
; summer starts May 1 [180 days]
; annual starts Jan 1 [365 days]
local nMsg, dimx, dimwy, ntim, save_FillValue, xts, nr, ns \
, nDay, nSea, iSea, mmdd, yyyy, mmddStrt, sdof, z1, smjo \
, spec, d, sm, pct, xtsVar, xvari
begin
; error check
nMsg = num( ismissing(x) )
if (nMsg.gt.0) then
print("mjo_spectra: currently, missing data not allowed: nMsg="+nMsg)
exit
end if
dimx = dimsizes(x)
dimwy = dimsizes(wy)
if (dimx(1).ne.dimwy(0)) then
print("band_pass_area_time: sizes of y/lat dimension do not match")
print(" dimsizes(wy)="+dimwy)
print(" dimx(1)=nlat="+dimx(1))
exit
end if
ntim = dimx(0)
xts = wgt_areaave_Wrap(x, wy, 1., 0) ; (time)
if (opt .and. isatt(opt,"detrend") .and. opt_at_detrend) then
if (isatt(xts,"_FillValue")) then
save_FillValue = x@_FillValue
delete(xts@_FillValue) ; avoid annoying warning msg
end if ; from dtrend
xts = dtrend(xts, False ) ; detrend
if (isvar("save_FillValue")) then
xts@_FillValue = save_FillValue ; reassign
end if
end if
nDay = 180
if (seaName.eq."winter") then
mmddStrt = 1101 ; Nov 1
end if
if (seaName.eq."summer") then
mmddStrt = 501 ; May 1
end if
if (seaName.eq."annual") then
nDay = 365
mmddStrt = 101 ; Jan 1
end if
yyyy = date/10000
mmdd = date - yyyy*10000
nSea = num(mmdd.eq.mmddStrt) ; number of seasons
iSea = ind(mmdd.eq.mmddStrt)
if (seaName.eq."winter") then
nSea = nSea-1 ; last season not complete
end if
;*****************************************************************
; calculate spectrum via segment averaging
; MJO Clivar says "no" to detrending/tapering
;*****************************************************************
if (opt .and. isatt(opt,"detrend_seg") ) then
d = opt_at_detrend_seg
else
d = 0 ; detrending opt: 0=>remove mean 1=>remove mean and detrend
end if
if (opt .and. isatt(opt,"smooth_seg") ) then
sm = opt_at_smooth_seg
else
sm = 0 ; smooth periodogram: 0 mean no smooth [pure periodgram]
end if
if (opt .and. isatt(opt,"taper_seg") ) then
pct = opt_at_taper_seg
else
pct = 0.0 ; percent tapered: (0.0 <= pct <= 1.0) 0.10 common.
end if
spec = new ( nDay/2, typeof(x))
spec = 0.0
z1 = 0.0
xtsVar= 0.0
do ns=0,nSea-1 ; segment averaging
iStrt = iSea(ns)
iLast = iSea(ns)+nDay-1
sdof = specx_anal(xts(iStrt:iLast),d,sm,pct)
spec = spec + sdof_at_spcx
xtsVar = xtsVar + sdof_at_xvari
z1 = z1 + 0.5*log((1+sdof_at_xlag1)/(1-sdof_at_xlag1)) ; Fischer z-transform
end do
spec = spec/nSea ; average spectra
z1 = z1/nSea ; mean z
xvari = xtsVar/nSea ; mean variance per segment
smjo = (/ sdof /)
smjo = 2*nSea ; # degrees freedom
smjo_at_bw = sdof_at_bw
smjo_at_spcx = spec
smjo_at_xvari = xvari
smjo_at_frq = sdof_at_frq
smjo_at_xlag1 = (exp(2*z1)-1)/(exp(2*z1)+1) ; match specx_anal
if (opt .and. isatt(opt,"smooth_seg") ) then
smjo = smjo*opt_at_smooth_seg
smjo_at_bw = (sdof_at_bw)*(opt_at_smooth_seg)
end if
return( smjo )
end
;*******************************************************************
undef("mjo_spectra")
procedure mjo_spectra \
(X[*][*][*]:numeric, date[*]:integer, wy[*]:numeric \
,latS[*]:numeric, latN[*]:numeric \
,lonL[*]:numeric, lonR[*]:numeric, nameRegion[*]:string \
,pltDir[1]:string, pltType[1]:string \
,pltName[1]:string, opt[1]:logical)
; driver for MJO CLIVAR segment averaging
; MJO CLIVAR: http://climate.snu.ac.kr/mjo_diagnostics/index.htm
; Plot spectra with x-axis as frequency with log scaling and
; y-axis as power times frequency
local frqmnPlt, frqmxPlt, nameSeason, nSeason, r, resP, tiXAxisString \
, ns, nr, plot, x, pltNameLocal, pltTypeLocal, pltPath, sMJO, splt\
, ifrqmxPlt, nRegion, ciLow, ciHigh
begin
frqmxPlt = 0.15 ; graphics X -axis [frequency]
frqmnPlt = 1./365.
nameSeason = (/ "winter", "summer", "annual" /)
nSeason = dimsizes(nameSeason)
nRegion = dimsizes(nameRegion)
;*******************************************************************
; general graphics resources
;*******************************************************************
r = True ; plot mods desired
r_at_gsnDraw = False ; do not draw
r_at_gsnFrame = False ; do not advance frame
r_at_tiYAxisString = "Power x freq" ; yaxis
r_at_xyLineThicknesses = (/2, 2, 1, 1/) ; Define line thicknesses
r_at_xyDashPatterns = (/0, 0, 1, 1/) ; Dash patterns
r_at_xyLineColors = (/"foreground","red","blue","blue"/)
r_at_vpHeightF = 0.4 ; change aspect ratio of plot
r_at_vpWidthF = 0.7
r_at_trXLog = True ; default
if (opt .and. isatt(opt,"logXAxis") .and. .not.opt_at_logXAxis) then
delete(r_at_trXLog) ; linear x (freq) axis
end if
if (opt .and. isatt(opt,"periodXAxis") .and. .not.opt_at_periodXAxis) then
tiXAxisString = "Frequency (cycles/day)" ; x (freq) axis
else
r_at_tmXBMode = "Explicit" ; default
r_at_tmXBValues = (/0.1, 0.05, 0.0333, 0.0167, 0.01, frqmnPlt/)
r_at_tmXBLabels = (/ 10, 20, 30 , 60 , 100, 365/)
tiXAxisString = "Period (days)" ; xaxis
end if
; linear Y is default
if (opt .and. isatt(opt,"logYAxis") .and. opt_at_logYAxis) then
r_at_trYLog = True
end if
resP = True
resP_at_gsnMaximize = True
resP_at_gsnStringFontHeightF = 0.025
resP_at_gsnPanelBottom = 0.05
plot = new ( nSeason, "graphic")
r_at_trXMaxF = frqmxPlt ; X Axis [rightmost freq]
if (pltType.eq."png") then
pltTypeLocal = "eps"
else
pltTypeLocal = pltType
end if
ciLow = 0.05 ; default
ciHigh = 0.95
if (opt .and. isatt(opt,"spcConLim")) then
if (dimsizes(opt_at_spcConLim).ne.2) then
print("mjo_spectra: spcConLim attribute must be size 2")
print(" dimsizes(spcConLim)="+dimsizes(opt_at_spcConLim))
exit
else
if (any(opt_at_spcConLim.lt.0 .or. opt_at_spcConLim.gt.1)) then
print("mjo_spectra: 0 < spcConLim < 1 violation")
print("spcConLim: "+opt_at_spcConLim)
exit
else
ciLow = opt_at_spcConLim(0) ; user specified bounds
ciHigh = opt_at_spcConLim(1)
end if
end if
end if
do nr=0,nRegion-1
; extract region
x = X(:,{latS(nr):latN(nr)},{lonL(nr):lonR(nr)})
wgty = wy({latS(nr):latN(nr)})
pltNameLocal = pltName+"_"+nameRegion(nr)
pltPath = pltDir+pltNameLocal
wks = gsn_open_wks(pltTypeLocal,pltPath)
do ns=0,nSeason-1
sMJO = mjo_spectra_season(x, date, wgty, nameSeason(ns), opt)
splt = specx_ci (sMJO, ciLow,ciHigh)
ifrqmxPlt= ind(sMJO_at_frq .ge. frqmnPlt .and. sMJO_at_frq .le. frqmxPlt)
r_at_gsnCenterString = nameSeason(ns)
if (isatt(r,"trXLog") .and. r_at_trXLog) then
r_at_trXMinF = frqmnPlt
end if
if (isatt(r,"trYLog") .and. r_at_trYLog) then
r_at_trYMinF = 0.95*min(splt(:,ifrqmxPlt))
r_at_trYMaxF = 1.05*max(splt(:,ifrqmxPlt))
end if
if (ns.eq.(nSeason-1)) then
r_at_tiXAxisString = tiXAxisString
end if
;plot(ns) = gsn_csm_xy(wks,sMJO_at_frq(ifrqmxPlt), splt(:,ifrqmxPlt),r)
work = splt
work = work*conform(work,sMJO_at_frq,1)
plot(ns) = gsn_csm_xy(wks,sMJO_at_frq(ifrqmxPlt), work(:,ifrqmxPlt),r)
if (ns.eq.(nSeason-1)) then
delete(r_at_tiXAxisString)
end if
delete(sMJO)
delete(splt)
delete(work)
delete(ifrqmxPlt)
end do
resP_at_txString = nameRegion(nr)+": "+x_at_long_name
gsn_panel(wks,plot,(/3,1/),resP)
if (pltType.eq."png") then
if (isatt(opt,"pltConvert")) then
pltConvert = opt_at_pltConvert ; convert options
else
pltConvert = " " ; default
end if
system("convert "+pltConvert+" "+pltPath+".eps "+pltPath+".png")
system("/bin/rm -f "+pltPath+".eps")
end if
delete(x)
delete(wgty)
end do
end
;*******************************************************************
undef("mjo_xcor_lag_season")
function mjo_xcor_lag_season(x[*]:numeric, y[*][*]:numeric, mxlag[1]:integer, opt[1]:logical)
; cross correlations at assorted lags
local dNam, dimy, N, yNew, xNew, y_Lead_x, x_Lead_y, ccr
begin
dNam = getvardims(y)
dimy = dimsizes(y)
N = dimy(1) ; N = mlon or nlat
yNew = y($dNam(1)$|:,$dNam(0)$|:)
xNew = conform(yNew,x,1)
y_Lead_x = esccr(xNew,yNew,mxlag) ; (N,mxlag+1)
x_Lead_y = esccr(yNew,xNew,mxlag)
ccr = new ( (/N,2*mxlag+1/), float)
ccr(:,0:mxlag-1) = x_Lead_y(:,1:mxlag:-1) ; "negative lag", -1 reverses order
ccr(:,mxlag:) = y_Lead_x(:,0:mxlag) ; "positive lag"
ccr = where(ccr.gt. 1.0, 1.0, ccr)
ccr = where(ccr.lt.-1.0,-1.0, ccr)
if (opt .and. isatt(opt,"smth9") ) then
if (abs(opt_at_smth9).eq.0.25) then
ccr = smth9(ccr, 0.50, opt_at_smth9, False)
else
print("mjo_xcor_lag_season: opt_at_smth9 must be 0.25 or -0.25")
print(" opt_at_smth9="+opt_at_smth9)
print(" NO SMOOTHING PERFORMED")
end if
end if
ccr!0 = dNam(1)
ccr&$dNam(1)$ = y&$dNam(1)$
ccr!1 = "lag"
ccr&lag = ispan(-mxlag,mxlag,1)
ccr&lag_at_units = "lag (days)"
ccr_at_long_name = "cross-correlation"
return(ccr(lag|:,$dNam(1)$|:))
end
;**********************************************************************
undef("mjo_xcor_lag_plot_ovly")
procedure mjo_xcor_lag_plot_ovly( \
ccr_a[*][*]:numeric, ccr_b[*][*]:numeric \
,pltType[1]:string, pltDir[1]:string \
,pltName[1]:string, opt[1]:logical)
local pltPath, pltTypeLocal, res1, res2, CCR_a, CCR_b, plot_a, plot_b
begin
pltPath = pltDir+"/"+pltName
if (pltType.eq."png") then
pltTypeLocal = "eps"
else
pltTypeLocal = pltType
end if
wks = gsn_open_wks(pltTypeLocal,pltPath)
if (opt .and. isatt(opt,"colorTable")) then
gsn_define_colormap(wks,opt_at_colorTable)
else
gsn_define_colormap(wks,"BlWhRe")
end if
res1 = True ; color precip
res1_at_gsnDraw = False
res1_at_gsnFrame = False
res1_at_gsnMaximize = True
res1_at_gsnPaperOrientation = "portrait"
res1_at_cnFillOn = True ; turn on color
res1_at_cnLinesOn = False
res1_at_gsnSpreadColors = True ; use full range of colormap
res1_at_cnLevelSelectionMode = "ManualLevels"; set manual contour levels
res1_at_cnMinLevelValF = -1.0 ; set min contour level
res1_at_cnMaxLevelValF = 1.0 ; set max contour level
res1_at_cnLevelSpacingF = 0.1 ; set contour spacing
res1_at_cnLabelBarEndLabelsOn= True
res1_at_cnLabelBarEndStyle = "ExcludeOuterBoxes"
res1_at_lbLabelAutoStride = True
res1_at_vpWidthF = 0.6 ; change aspect ratio of plot
res1_at_vpHeightF = 0.4
;res1_at_gsnMajorLatSpacing = 10 ; change maj lat tm spacing
res1_at_tiYAxisString = "lag (days)"
if (isatt(opt,"tiMainString")) then
res1_at_tiMainString = opt_at_tiMainString
end if
if (isatt(opt,"gsnLeftString")) then
res1_at_gsnLeftString = opt_at_gsnLeftString
end if
if (isatt(opt,"gsnCenterString")) then
res1_at_gsnCenterString = opt_at_gsnCenterString
end if
if (isatt(opt,"gsnRightString")) then
res1_at_gsnRightString = opt_at_gsnRightString
end if
;************************************************
; resource list for second data array
;************************************************
res2 = True ; U
res2_at_gsnDraw = False
res2_at_gsnFrame = False
res2_at_cnLevelSelectionMode = "ManualLevels" ; set manual contour levels
res2_at_cnMinLevelValF = -1.0
res2_at_cnMaxLevelValF = 1.0
res2_at_cnLevelSpacingF = 0.1
res2_at_cnLineLabelsOn = True
res2_at_gsnContourZeroLineThicknessF = 0. ; Eliminate 0 line
res2_at_gsnContourNegLineDashPattern = 1 ; negative contours dash pattern
res2_at_cnLineThicknessF = 2. ; line thickness
res2_at_cnInfoLabelOn = False
CCR_a = ccr_a ; possible smooth and delete of attribute
CCR_b = ccr_b
if (opt .and. isatt(opt,"smth9") .and. abs(opt_at_smth9).eq.0.25) then
CCR_a = smth9(CCR_a, 0.50, opt_at_smth9, False)
end if
delete(CCR_a_at_long_name)
plot_a = gsn_csm_contour(wks,CCR_a,res1) ; contour the variable
if (opt .and. isatt(opt,"smth9") .and. abs(opt_at_smth9).eq.0.25) then
CCR_b = smth9(CCR_b, 0.50, opt_at_smth9, False)
end if
delete(CCR_b_at_long_name)
plot_b = gsn_csm_contour(wks,CCR_b,res2) ; contour the variable
overlay(plot_a, plot_b)
draw(plot_a)
frame(wks)
if (pltType.eq."png") then
if (isatt(opt,"pltConvert")) then
pltConvert = opt_at_pltConvert ; convert options
else
pltConvert = " " ; default
end if
system("convert "+pltConvert+" "+pltPath+".eps "+pltPath+".png")
system("/bin/rm -f "+pltPath+".eps")
end if
end
;**********************************************************************
undef("mjo_xcor_lag_ovly_panel")
procedure mjo_xcor_lag_ovly_panel( \
ccr_a[*][*]:numeric, ccr_b[*][*]:numeric \
,ccr_c[*][*]:numeric, ccr_d[*][*]:numeric \
,pltType[1]:string, pltDir[1]:string \
,pltName[1]:string, opt[1]:logical)
local pltPath, pltTypeLocal, res1, res2, CCR1, CCR2 \
, plt1, plt2, resP, plot
begin
pltPath = pltDir+"/"+pltName
if (pltType.eq."png") then
pltTypeLocal = "eps"
else
pltTypeLocal = pltType
end if
wks = gsn_open_wks(pltTypeLocal,pltPath)
if (opt .and. isatt(opt,"colorTable")) then
gsn_define_colormap(wks,opt_at_colorTable)
else
gsn_define_colormap(wks,"BlWhRe")
end if
res1 = True ; color precip
res1_at_gsnDraw = False
res1_at_gsnFrame = False
res1_at_gsnMaximize = True
res1_at_gsnPaperOrientation = "portrait"
res1_at_cnFillOn = True ; turn on color
res1_at_cnLinesOn = False
res1_at_gsnSpreadColors = True ; use full range of colormap
res1_at_cnLevelSelectionMode = "ManualLevels"; set manual contour levels
res1_at_cnMinLevelValF = -1.0 ; set min contour level
res1_at_cnMaxLevelValF = 1.0 ; set max contour level
res1_at_cnLevelSpacingF = 0.1 ; set contour spacing
res1_at_cnLabelBarEndLabelsOn= True
res1_at_cnLabelBarEndStyle = "ExcludeOuterBoxes"
res1_at_cnInfoLabelOn = False
res1_at_lbLabelBarOn = False ; turn off individual cb's
res1_at_vpWidthF = 0.6 ; change aspect ratio of plot
res1_at_vpHeightF = 0.4
;res1_at_gsnMajorLatSpacing = 10 ; change maj lat tm spacing
res1_at_tiYAxisString = "lag (days)"
if (isatt(opt,"tiMainString")) then
res1_at_tiMainString = opt_at_tiMainString
end if
if (isatt(opt,"gsnLeftString")) then
res1_at_gsnLeftString = opt_at_gsnLeftString
end if
if (isatt(opt,"gsnCenterString")) then
res1_at_gsnCenterString = opt_at_gsnCenterString
end if
if (isatt(opt,"gsnRightString")) then
res1_at_gsnRightString = opt_at_gsnRightString
end if
;************************************************
; resource list for second data array
;************************************************
res2 = True ; U
res2_at_gsnDraw = False
res2_at_gsnFrame = False
res2_at_cnLevelSelectionMode = "ManualLevels" ; set manual contour levels
res2_at_cnMinLevelValF = -1.0
res2_at_cnMaxLevelValF = 1.0
res2_at_cnLevelSpacingF = 0.1
res2_at_cnLineLabelsOn = True
res2_at_gsnContourZeroLineThicknessF = 0. ; Eliminate 0 line
res2_at_gsnContourNegLineDashPattern = 1 ; negative contours dash pattern
;res2_at_cnLineThicknessF = 2. ; line thickness
res2_at_cnInfoLabelOn = False
plot = new( 2, "graphic")
do np=0,1
if (np.eq.0) then
CCR1 = ccr_a ; possible smooth and delete of attribute
CCR2 = ccr_b
else
CCR1 = ccr_c ; possible smooth and delete of attribute
CCR2 = ccr_d
end if
if (opt .and. isatt(opt,"smth9") .and. abs(opt_at_smth9).eq.0.25) then
CCR1 = smth9(CCR1, 0.50, opt_at_smth9, False)
CCR2 = smth9(CCR2, 0.50, opt_at_smth9, False)
end if
delete(CCR1_at_long_name)
plt1 = gsn_csm_contour(wks,CCR1,res1) ; contour the variable
delete(CCR2_at_long_name)
plt2 = gsn_csm_contour(wks,CCR2,res2) ; contour the variable
overlay(plt1, plt2)
plot(np) = plt1
delete(CCR1) ; size may change
delete(CCR2)
end do
resP = True ; panel
resP_at_lbLabelAutoStride = True
resP_at_gsnPanelLabelBar = True ; add common colorbar
;resP_at_lbLabelFontHeightF = 0.007 ; make labels smaller
if (isatt(opt,"txString")) then
resP_at_txString = opt_at_txString
end if
gsn_panel(wks,plot,(/2,1/),resP) ; now draw as one plot
if (pltType.eq."png") then
if (isatt(opt,"pltConvert")) then
pltConvert = opt_at_pltConvert ; convert options
else
pltConvert = " " ; default
end if
system("convert "+pltConvert+" "+pltPath+".eps "+pltPath+".png")
system("/bin/rm -f "+pltPath+".eps")
end if
end
;***********************************************************
undef("mjo_xcor_lag")
function mjo_xcor_lag (xio[*]:numeric , y[*][*]:numeric \
,date[*]:integer, mxlag[1]:integer, seaName:string, opt[1]:logical)
; driver to calculate mean seasonal lags
local ns, ntim, dimy, ntmy, my, nDay, mmddStrt, yyyy, mmdd, nSea, iSea, rtyp
, rLag, nLag, rr, iStrt, iLast
begin
ntim = dimsizes(xio)
dimy = dimsizes(y)
ntmy = dimy(0)
my = dimy(1)
if (ntim.ne.ntmy) then
print("mjo_xcor_lag: time dimensions do not match")
print(" ntim="+ntim+" ntmy="+ntmy)
exit
end if
nDay = 180
if (seaName.eq."winter") then
mmddStrt = 1101 ; Nov 1
end if
if (seaName.eq."summer") then
mmddStrt = 501 ; May 1
end if
if (seaName.eq."annual") then
nDay = 365
mmddStrt = 101 ; Jan 1
end if
yyyy = date/10000
mmdd = date - yyyy*10000
nSea = num(mmdd.eq.mmddStrt) ; number of 'seaNam' seasons
iSea = ind(mmdd.eq.mmddStrt)
if (seaName.eq."winter") then
nSea = nSea-1 ; last season not complete
end if
rtyp = "float"
if (typeof(xio).eq."double" .or. typeof(y).eq."double") then
rtyp = "double"
end if
nLag = 2*mxlag+1
rz = new ( (/nLag,my/), rtyp, getFillValue(y) )
rz = 0.0
do ns=0,nSea-1 ; seasonal averaging
iStrt = iSea(ns)
iLast = iSea(ns)+nDay-1
if (iLast.gt.(ntim-1)) then
break
end if
rLag = mjo_xcor_lag_season (xio(iStrt:iLast), y(iStrt:iLast,:) ,mxlag, opt)
rz = rz + 0.5*log((1+rLag)/(1-rLag)) ; Fischer z-transform
end do
rz = rz/nSea ; mean Fischer-z
rz = (exp(2*rz)-1)/(exp(2*rz)+1) ; transform back
copy_VarMeta(rLag, rz)
return( rz )
end
; ====================================================
undef("mjo_wavenum_freq_season")
function mjo_wavenum_freq_season \
(X[*][*]:numeric, date[*]:integer \
,seaName[1]:string, opt:logical)
;
; For each 'seaName' (say, 'winter') over a number
; of different years, calculate the
; pooled space-time spectra
; MJO CLIVAR: wavenumber-frequency spectra
; winter starts Nov 1 [180 days]
; summer starts May 1 [180 days]
; annual starts Jan 1 [365 days]
local nMsg, dimX, ntim, mlon, x, N, pow, xAvg, xSeason \
, nDay, ny, nYear, iSea, mmdd, yyyy, mmddStrt, work, work1\
, d, sm, pct, xVarSea, kSea, wave, freq, iStrt, iLast
begin
; error check
nMsg = num( ismissing(X) )
if (nMsg.gt.0) then
print("mjo_wavenum_freq_season: currently, missing data not allowed: nMsg="+nMsg)
exit
end if
dimX = dimsizes(X)
ntim = dimX(0) ; total times
mlon = dimX(1)
x = X ; local copy (time,lon)
; x may change
if (isatt(x,"_FillValue")) then
delete(x@_FillValue) ; avoid annoying warning msg
end if ; from dtrend
nDay = 180
if (seaName.eq."winter") then
mmddStrt = 1101 ; Nov 1
end if
if (seaName.eq."summer") then
mmddStrt = 501 ; May 1
end if
if (seaName.eq."annual") then
nDay = 365
mmddStrt = 101 ; Jan 1
end if
yyyy = date/10000
mmdd = date - yyyy*10000
nYear = num(mmdd.eq.mmddStrt) ; number of seasons
iSea = ind(mmdd.eq.mmddStrt)
;*****************************************************************
; For a specific season, calculate spectra via averaging
; over each seasonal segment.
; MJO Clivar says "no" to detrending/tapering.
; Hence, the following are just 'place holders'
;*****************************************************************
x = dtrend_leftdim(x, False ) ; detrend overall series in time
work = new((/2,mlon ,nDay /), typeof(x), "No_FillValue")
; the +1 is for the 0-th component
pow = new((/ mlon+1,nDay+1/), typeof(x), "No_FillValue")
pow!0 = "wave"
pow!1 = "frq"
pow = 0.0
xAvgSea = 0.0
xVarSea = 0.0 ; variance (raw)
xVarTap = 0.0 ; variance after tapering
kSea = 0 ; count os seasons used
N = nDay ; convenience
if (opt .and. isatt(opt,"taper_seg") ) then
pct = opt_at_taper_seg ; over-ride default
else
pct = 0.10 ; default; taper 10% of series
end if
; -------------------------------------------------------------------
; NCL does not have a complex 2D FFT at this time.
; Perform a "poorman's" complex 2D FFT by looping over time and space.
; Also, makes the code more clear (IMHO)
; -------------------------------------------------------------------
do ny=0,nYear-1 ; season (segment) averaging
iStrt = iSea(ny) ; start index for current season
iLast = iSea(ny)+nDay-1 ; last
if (iLast.gt.(ntim-1)) then
break
end if
;print("iStrt="+iStrt+" ("+yyyy(iStrt)+mmdd(iStrt)+"); " \
; +"iLast="+iLast+" ("+yyyy(iLast)+mmdd(iLast)+")" )
xSeason = x(iStrt:iLast,:) ; keep meta data
xAvg = avg( xSeason ) ; season average all time/lon
xSeason = xSeason - xAvg ; remove season time-lon mean
xVarSea = xVarSea + variance( xSeason ) ; overall variance
kSea = kSea + 1
xSeason = taper_leftdim( xSeason, pct, 0) ; for each lon detrend in time
xVarTap = xVarTap + variance( xSeason ) ; variance after tapering
do nt=0,nDay-1 ; each time perform complex fft(lon)
work(:,:,nt) = cfftf( xSeason(nt,:), 0.0, 0) ; space
end do
work = work/mlon ; normalize by # lon samples
do ml=0,mlon-1 ; each lon perform complex fft(time)
work(:,ml,:) = cfftf( work(0,ml,:), work(1,ml,:), 0) ; time
end do
work = work/nDay ; normalize by # time samples
; work(2,wave,freq)
; unscramble fft+ calculate power
pow = pow + resolveWavesHayashi( work, nDay, 1) ; (wave,freq)
end do
xVarSea = xVarSea/kSea ; pooled seasonal variance
xVarTap = xVarTap/kSea ; pooled seasonal variance
pow = pow/kSea ; pooled spectra
if (isatt(X,"long_name")) then
pow_at_long_name = x_at_long_name
end if
; add meta data for use upon return
pow!0 = "wave"
pow!1 = "freq"
wave = ispan(-mlon/2,mlon/2,1)
wave!0 = "wave"
wave_at_long_name= "zonal wavenumber" ; for plot
wave&wave = wave
freq = fspan(-1*nDay/2,nDay/2,nDay+1)/nDay
freq!0 = "freq"
freq_at_long_name= "frequency (cycles/day)" ; for plot
freq_at_units= "cycles/day"
freq&freq = freq
pow&wave = wave
pow&freq = freq
return( pow )
end
;--------------------------------------------------------------------
undef("mjo_wavenum_freq_season_plot")
procedure mjo_wavenum_freq_season_plot (wf[*][*],seaName[1]:string\
,pltDir[1]:string, pltType[1]:string \
,pltName[1]:string, opt[1]:logical)
begin
pltPath= pltDir+"/"+pltName+"."+seaName
if (pltType.eq."png") then
pltTypeLocal = "eps"
else
pltTypeLocal = pltType
end if
wks = gsn_open_wks(pltTypeLocal, pltPath)
if (opt .and. isatt(opt,"colorTable")) then
gsn_define_colormap(wks,opt_at_colorTable)
else
gsn_define_colormap(wks,"prcp_2")
end if
res = True ; plot mods desired
res_at_gsnFrame = False
;res_at_gsnDraw = False
res_at_cnFillOn = True ; turn on color
res_at_gsnSpreadColors = True ; use full range of colormap
res_at_lbLabelAutoStride = True
if (isatt(opt,"tiMainString")) then
res_at_tiMainString = opt_at_tiMainString
else
res_at_tiMainString = changeCase(seaName, "up")
end if
if (isatt(opt,"gsnLeftString")) then
res_at_gsnLeftString = opt_at_gsnLeftString
end if
if (isatt(opt,"gsnCenterString")) then
res_at_gsnCenterString = opt_at_gsnCenterString
end if
if (isatt(opt,"gsnRightString")) then
res_at_gsnRightString = opt_at_gsnRightString
end if
if (isatt(opt,"cnLinesOn")) then
res_at_cnLinesOn = opt_at_cnLinesOn
if (.not.res_at_cnLinesOn) then
res_at_cnLineLabelsOn = False
end if
end if
if (isatt(opt,"cnLevelSelectionMode")) then
res_at_cnLevelSelectionMode = opt_at_cnLevelSelectionMode
end if
if (isatt(opt,"cnLevelSelectionMode")) then
res_at_cnLevelSelectionMode = opt_at_cnLevelSelectionMode
end if
if (isatt(opt,"cnMinLevelValF")) then
res_at_cnMinLevelValF = opt_at_cnMinLevelValF
end if
if (isatt(opt,"cnMaxLevelValF")) then
res_at_cnMaxLevelValF = opt_at_cnMaxLevelValF
end if
if (isatt(opt,"cnLevelSpacingF")) then
res_at_cnLevelSpacingF = opt_at_cnLevelSpacingF
end if
if (isatt(opt,"cnLevels")) then
res_at_cnLevels = opt_at_cnLevels
end if
NW = 6
if (isatt(opt,"maxWavePlot")) then ; wave [Y] axis
NW = opt_at_maxWavePlot
end if
fMin = -0.05
fMax = 0.05
if (isatt(opt,"minFreqPlot")) then
fMin = opt_at_minFreqPlot
end if
if (isatt(opt,"maxFreqPlot")) then
fMax = opt_at_maxFreqPlot
end if
;res_at_trXMinF = fMin
;res_at_trXMaxF = fMax
WORK = wf({0:NW},{fMin:fMax})
if (opt .and. isatt(opt,"smth9") .and. opt_at_smth9) then
WORK = smth9(WORK, 0.50, 0.25, False)
end if
izero = ind(WORK&freq .eq. 0.0)
WORK(:,izero) = min(WORK) ; 0th freq
plot = gsn_csm_contour(wks, WORK , res)
resp = True ; polyline mods desired
;resp_at_gsLineThicknessF = 2.0 ; thickness of lines
resp_at_gsLineDashPattern= 11
tres = True
tres_at_txFontHeightF = 0.0175
day = 30
fline = 1./day
;resp_at_gsLineLabelString= day+"d" ; adds a line label string
gsn_polyline(wks,plot,(/fline,fline/),(/ 0.,NW/), resp)
gsn_text(wks,plot, (day+"d"),fline+0.005,0.93*NW,tres)
day = 80
fline = 1/day
;resp_at_gsLineLabelString= day+"d" ; adds a line label string
gsn_polyline(wks,plot,(/fline,fline/),(/ 0.,NW/), resp)
gsn_text(wks,plot, (day+"d"),fline+0.005,0.93*NW,tres)
frame(wks)
if (pltType.eq."png") then
if (isatt(opt,"pltConvert")) then
pltConvert = opt_at_pltConvert ; convert options
else
pltConvert = " " ; default
end if
system("convert "+pltConvert+" "+pltPath+".eps "+pltPath+".png")
system("/bin/rm -f "+pltPath+".eps")
end if
if (pltType.eq."x11" .and. isatt(opt,"debug") .and. opt_at_debug) then
res_at_gsnCenterString = "wf"
plot = gsn_csm_contour(wks, wf, res)
res_at_gsnCenterString = "wf({-5:5},{-0.075:0.075})"
plot = gsn_csm_contour(wks, wf({-5:5},{-0.075:0.075}), res)
end if
end
;----------------------------------------------------------------
undef("mjo_cross")
function mjo_cross (X[*][*][*]:numeric, Y[*][*][*]:numeric \
,segLen[1]:integer , segOverLap[1]:integer \
,opt:logical)
local nMsgX, nMsgY, dimX, dimY, ntim, nlat, mlon, switc \
, x, y, pct, STC, stc, namStc, kseg, freq, wave \
, ntStrt, ntLast
begin
; error check
dimX = dimsizes(X)
dimY = dimsizes(Y)
if (.not.all(dimX.eq.dimY)) then
print("mjo_cross: X and Y must be same size")
print(" dimX="+dimX+" dimY="+dimY)
exit
end if
nMsgX = num( ismissing(X) )
nMsgY = num( ismissing(Y) )
;if ((nMsgX+nMsgY).gt.0) then ; checked in builin function
; print("mjo_cross: currently, missing data not allowed: nMsgX="+nMsgX+" nMsgY="+nMsgY)
; exit
;end if
ntim = dimX(0) ; total times
nlat = dimX(1)
mlon = dimX(2)
x = X ; local copy (time,lat,lon)
y = Y
x = decompose2SymAsym( x, 0) ; decompose
y = decompose2SymAsym( y, 0)
; x may change
if (isatt(x,"_FillValue")) then
delete(x@_FillValue) ; avoid annoying warning msg
end if ; from dtrend
if (isatt(y,"_FillValue")) then
delete(y@_FillValue)
end if
x = dtrend_leftdim(x, False ) ; detrend overall series in time
y = dtrend_leftdim(y, False )
if (opt .and. isatt(opt,"taper_seg") ) then
pct = opt_at_taper_seg ; over-ride default
else
pct = 0.10 ; default; taper 10% of series
end if
x = taper_leftdim(x, pct, 0 ) ; taper in time
y = taper_leftdim(y, pct, 0 )
; -------------------------------------------------------------------
STC = new ( (/16,segLen/2+1,mlon+1/), "double", 1d20)
STC = 0.0
kseg = 0
ntStrt = 0
switch = True
do while (switch)
ntLast = ntStrt + segLen-1
if (ntLast.gt.(ntim-1)) then
switch = False
break
end if
XX = taper_leftdim( x(ntStrt:ntLast,:,:), pct, 0)
YY = taper_leftdim( y(ntStrt:ntLast,:,:), pct, 0)
kseg = kseg+1
STC = STC + mjo_cross_segment(XX,YY,0)
ntStrt = ntStrt + abs(segOverLap)
end do
STC = STC/kseg ; average
; use averaged power
mjo_cross_coh2pha(STC, 0) ; to get coh2 / phase
freq = fspan(0,0.5,segLen/2+1)
freq_at_long_term = " frq (cycles per day)" ; for plot
freq_at_units = "cycles per day"
freq!0 = "freq"
wave = fspan(-mlon/2,mlon/2,mlon+1)
wave_at_long_term = "zonal wavenumber"
wave_at_units = "zonal wavenumber"
wave!0 = "wave"
STC!0 = "var"
STC!1 = "freq"
STC!2 = "wave"
STC&freq = freq
STC&wave = wave
STC_at_varName = (/"PEE1_SYM" , "PEE1_ASY" \
,"PEE2_SYM" , "PEE2_ASY" \
,"P12_SYM" , "P12_ASY" \
,"Q12_SYM" , "Q12_ASY" \
,"COH2_SYM" , "COH2_ASY" \
,"PHAS_SYM" , "PHAS_ASY" \
,"V1_SYM" , "V1_ASY" \
,"V2_SYM" , "V2_ASY" /)
STC_at_segmentLength = segLen
STC_at_segmentOverLap = segOverLap
STC_at_segmentRepeat = segLen-segOverLap
STC_at_number_of_segments = kseg
STC_at_dof = 2.667*kseg ; conservative estimate
; 2.667 is for 1-2-1 smoother
p = (/0.80, 0.85, 0.90, 0.925, 0.95, 0.99/); probability levels
STC_at_prob = p
STC_at_prob_coh2 = 1.-(1.-p^(0.5*STC_at_dof -1))
return( STC )
end
; -----------
undef("addHorVertLinesCross")
function addHorVertLinesCross(wks[1]:graphic, plot[1]:graphic \
,nw[1], dumy[*]:graphic )
; freq [y] axis: Add horizontal lines that explicitly
; print time in days. This assumes the units
; of the freq axis are "cpd" [cycles per day]
local gsres, txres, xx, dely, m\nwl, nwr
begin
gsres = True
gsres_at_gsLineDashPattern = 1
nwl = -nw+3.5 ; left
nwr = nw ; right
dumy(0) = gsn_add_polyline(wks, plot, (/ 0, 0/),(/ 0.0 , 0.5 /),gsres)
dumy(1) = gsn_add_polyline(wks, plot, (/nwl,nwr/),(/ 1./80, 1./80/),gsres)
dumy(2) = gsn_add_polyline(wks, plot, (/nwl,nwr/),(/ 1./20, 1./20/),gsres)
dumy(3) = gsn_add_polyline(wks, plot, (/nwl,nwr/),(/ 1./10, 1./10/),gsres)
dumy(4) = gsn_add_polyline(wks, plot, (/nwl,nwr/),(/ 1./5 , 1./5 /),gsres)
dumy(5) = gsn_add_polyline(wks, plot, (/nwl,nwr/),(/ 1./3 , 1./3 /),gsres)
txres = True
txres_at_txJust = "CenterLeft"
txres_at_txFontHeightF = 0.013
xx = -nw+0.3
dely = 0.000 ; yy
dumy(6) = gsn_add_text(wks, plot, "3 days" , xx ,(1./3 +dely),txres)
dumy(7) = gsn_add_text(wks, plot, "5 days" , xx ,(1./5 +dely),txres)
dumy(8) = gsn_add_text(wks, plot, "10 days", xx ,(1./10+dely),txres)
dumy(9) = gsn_add_text(wks, plot, "20 days", xx ,(1./20+dely),txres)
dumy(10)= gsn_add_text(wks, plot, "80 days", xx ,(1./80+dely),txres)
return(plot)
end
; -----------
undef("mjo_elimIsolatedValues")
procedure mjo_elimIsolatedValues(c2, ph1, ph2, iopt)
; cosmetic and arbitrary ... eliminate isolated points ...
; crude and not complete but 'good enough'
local npts, dimc2, nfreq, nwave, nf, nw
begin
i = iopt
if (iopt.le.0) then
i = 1
end if
npts = 0
dimc2 = dimsizes(c2)
nfreq = dimc2(0)
nwave = dimc2(1)
do nf=i,nfreq-i-1
do nw=i,nwave-i-1
if (.not.ismissing(c2(nf,nw))) then
if (all(ismissing(c2(nf-i,nw )).and.ismissing(c2(nf+i,nw )) .and.\
ismissing(c2(nf ,nw-i)).and.ismissing(c2(nf ,nw+i)))) then
npts = npts + 1
c2(nf,nw) = c2@_FillValue
ph1(nf,nw) = ph1@_FillValue
ph2(nf,nw) = ph2@_FillValue
end if
end if
end do
nw = 0 ; left (freq) boundary
if (.not.ismissing(c2(nf,nw))) then
if (all(ismissing(c2(nf-i,nw )).and.ismissing(c2(nf+i,nw )) .and. \
ismissing(c2(nf ,nw+1)))) then
npts = npts + 1
c2(nf,nw) = c2@_FillValue
ph1(nf,nw) = ph1@_FillValue
ph2(nf,nw) = ph2@_FillValue
end if
end if
nw = nwave-1 ; right (freq) boundary
if (.not.ismissing(c2(nf,nw))) then
if (all(ismissing(c2(nf-1,nw )).and.ismissing(c2(nf+i,nw )) .and. \
ismissing(c2(nf ,nw-1)))) then
npts = npts + 1
c2(nf,nw) = c2@_FillValue
ph1(nf,nw) = ph1@_FillValue
ph2(nf,nw) = ph2@_FillValue
end if
end if
end do
;print("***** elimIsolatedValues: npts = "+npts+" set to _FillValue")
end
; -----------
undef("mjo_cross_plot")
procedure mjo_cross_plot(STC[*][*][*]:numeric \
,pltDir[1]:string, pltType[1]:string \
,pltName[1]:string, opt[1]:logical )
local pltPath, pltTypeLocal,wks, res
begin
pltPath= pltDir+"/"+pltName
if (pltType.eq."png") then
pltTypeLocal = "eps"
else
pltTypeLocal = pltType
end if
wks = gsn_open_wks(pltTypeLocal, pltPath)
if (opt .and. isatt(opt,"colorTable")) then
gsn_define_colormap(wks,opt_at_colorTable)
else
gsn_define_colormap(wks,"amwg")
end if
plot = new ( 2, "graphic")
res = True ; plot mods desired
res_at_gsnDraw = False
res_at_gsnFrame = False
res_at_cnFillOn = True ; turn on color
res_at_cnFillMode = "RasterFill" ; match WMO Clivar
res_at_gsnSpreadColors = True ; use full range of colormap
if (opt .and. isatt(opt,"gsnSpreadColorStart") ) then
res_at_gsnSpreadColorStart = opt_at_gsnSpreadColorStart
end if
if (opt .and. isatt(opt,"gsnSpreadColorEnd") ) then
res_at_gsnSpreadColorEnd = opt_at_gsnSpreadColorEnd
end if
res_at_cnLinesOn = False
res_at_cnLineLabelsOn = False
res_at_cnLevelSelectionMode = "ManualLevels"
res_at_cnMinLevelValF = 0.05
res_at_cnMaxLevelValF = 0.65 ; correlation^2 = 0.8
res_at_cnLevelSpacingF = 0.05
res_at_cnInfoLabelOn = False
res_at_lbLabelBarOn = False ; no individual label bars
if (.not.opt .or. .not.isatt(opt,"pltPhase") .or. opt_at_pltPhase) then
plotPhase = True
res_at_vcRefMagnitudeF = 1.0 ; define vector ref mag
res_at_vcRefLengthF = 0.01 ; define length of vec ref
res_at_vcRefAnnoOrthogonalPosF = -1.0 ; move ref vector
res_at_vcRefAnnoArrowLineColor = "black" ; change ref vector color
res_at_vcMinDistanceF = 0.0075 ; thin out vectors
res_at_vcMapDirection = False ;
res_at_vcRefAnnoOn = False ; do not draw
res_at_gsnScalarContour = True ; contours desired
else
plotPhase = False
end if
res_at_gsnLeftString = "Coh^2: Symmetric"
res_at_gsnRightString = "10%="+sprintf("%3.2f", STC_at_prob_coh2(2))+" "\
+" 5%="+sprintf("%3.2f", STC_at_prob_coh2(4))
if (opt .and. isatt(opt,"pltZonalWaveNumber") ) then
nWavePlt = opt_at_pltZonalWaveNumber
else
nWavePlt = 15 ; default
end if
;---------------------------------------------------------------
; dispersion: curves
;---------------------------------------------------------------
rlat = 0.0
Ahe = (/50.,25.,12./)
nWaveType = 6
nPlanetaryWave = 50
nEquivDepth = dimsizes(Ahe)
Apzwn = new((/nWaveType,nEquivDepth,nPlanetaryWave/),"double",1e20)
Afreq = Apzwn
genDispersionCurves(nWaveType, nEquivDepth, nPlanetaryWave, rlat, Ahe, Afreq, Apzwn )
;---------------------------------------------------------------
; dispersion curve and text plot resources
;---------------------------------------------------------------
dcres = True
dcres_at_gsLineThicknessF = 2.0
dcres_at_gsLineDashPattern = 0
txres = True
;txres_at_txJust = "CenterLeft"
txres_at_txPerimOn = True
txres_at_txFontHeightF = 0.013
txres_at_txBackgroundFillColor = "Background"
;---------------------------------------------------------------
; plot symmetric data
;---------------------------------------------------------------
n = 8
c2s = STC(n,:,{-nWavePlt:nWavePlt})
c2s@_FillValue = 1e20
c2s(0,:) = c2s@_FillValue
c2s = where(c2s.lt.0.05, c2s@_FillValue, c2s) ; mask
n = 12
phs1 = STC(n,:,{-nWavePlt:nWavePlt})
phs1_at_long_name = "symmetric phase-1"
phs1@_FillValue = c2s@_FillValue
phs1(0,:) = phs1@_FillValue
phs1 = where(c2s.lt.0.05, phs1@_FillValue, phs1) ; mask
n = 14
phs2 = STC(n,:,{-nWavePlt:nWavePlt})
phs2_at_long_name = "symmetric phase-2"
phs2@_FillValue = c2s@_FillValue
phs2(0,:) = phs2@_FillValue
phs2 = where(c2s.lt.0.05, phs2@_FillValue, phs2) ; mask
if (opt .and. isatt(opt,"pltProb") ) then
np = ind(c2s_at_prob_coh2 .eq. opt_at_pltProb)
if (.not.ismissing(np)) then
c2s = where(c2s.lt.STC_at_prob_coh2(np), c2s@_FillValue, c2s)
phs1 = where(ismissing(c2s), phs1@_FillValue, phs1)
phs2 = where(ismissing(c2s), phs2@_FillValue, phs2)
end if
end if
if (opt .and. isatt(opt,"coh2Cutoff") ) then
c2s = where(c2s.lt.opt_at_coh2Cutoff, c2s@_FillValue, c2s)
end if
if (opt .and. isatt(opt,"phaseCutoff") ) then
phs1 = where(c2s.lt.opt_at_phaseCutoff, phs1@_FillValue, phs1)
phs2 = where(c2s.lt.opt_at_phaseCutoff, phs2@_FillValue, phs2)
end if
if (opt .and. isatt(opt,"elimIsoVals") .and. .not.opt_at_elimIsoVals) then
print("mjo_cross_plot: no values eliminated")
else
mjo_elimIsolatedValues (c2s, phs1, phs2, 1)
mjo_elimIsolatedValues (c2s, phs1, phs2, 2)
end if
;res_at_gsnCenterString = "var="+n+" "+STC_at_varName(n)
if (plotPhase) then
scl_one = sqrt(1./(phs1^2 + phs2^2))
phs1 = scl_one*phs1
phs2 = scl_one*phs2
plot(0) = gsn_csm_vector_scalar(wks, phs1, phs2, c2s, res)
else
plot(0) = gsn_csm_contour(wks, c2s, res)
end if
dums = new (25, "graphic") ; for _add_ *s*ymmetric graphical objects
plot(0) = addHorVertLinesCross(wks, plot(0), nWavePlt, dums)
dumdcs = new ( 25, "graphic") ; dispersion curves symmetric (dcs)
dumdcs( 0) = gsn_add_polyline(wks,plot(0),Apzwn(3,0,:),Afreq(3,0,:),dcres)
dumdcs( 1) = gsn_add_polyline(wks,plot(0),Apzwn(3,1,:),Afreq(3,1,:),dcres)
dumdcs( 2) = gsn_add_polyline(wks,plot(0),Apzwn(3,2,:),Afreq(3,2,:),dcres)
dumdcs( 3) = gsn_add_polyline(wks,plot(0),Apzwn(4,0,:),Afreq(4,0,:),dcres)
dumdcs( 4) = gsn_add_polyline(wks,plot(0),Apzwn(4,1,:),Afreq(4,1,:),dcres)
dumdcs( 5) = gsn_add_polyline(wks,plot(0),Apzwn(4,2,:),Afreq(4,2,:),dcres)
dumdcs( 6) = gsn_add_polyline(wks,plot(0),Apzwn(5,0,:),Afreq(5,0,:),dcres)
dumdcs( 7) = gsn_add_polyline(wks,plot(0),Apzwn(5,1,:),Afreq(5,1,:),dcres)
dumdcs( 8) = gsn_add_polyline(wks,plot(0),Apzwn(5,2,:),Afreq(5,2,:),dcres)
dumdcs( 9) = gsn_add_text(wks,plot(0),"Kelvin",11.5,.40,txres)
dumdcs(10) = gsn_add_text(wks,plot(0),"n=1 ER",-10.7,.07,txres)
dumdcs(11) = gsn_add_text(wks,plot(0),"n=1 IG",-3.0,.45,txres)
dumdcs(12) = gsn_add_text(wks,plot(0),"h=50",-14.0,.78,txres)
dumdcs(13) = gsn_add_text(wks,plot(0),"h=25",-14.0,.60,txres)
dumdcs(14) = gsn_add_text(wks,plot(0),"h=12",-14.0,.46,txres)
;---------------------------------------------------------------
; plot asymmetric data
;---------------------------------------------------------------
n = 9
res_at_gsnLeftString = "Coh^2: Asymmetric"
res_at_gsnRightString = "10%="+sprintf("%3.2f", STC_at_prob_coh2(2))+" "\
+" 5%="+sprintf("%3.2f", STC_at_prob_coh2(4))
c2a = STC(n,:,{-nWavePlt:nWavePlt})
c2a@_FillValue = 1e20
c2a(0,:) = c2a@_FillValue
c2a = where(c2a.lt.0.05, c2a@_FillValue, c2a) ; mask
n = 13
pha1 = STC(n,:,{-nWavePlt:nWavePlt})
pha1_at_long_name = "asymmetric phase-1"
pha1@_FillValue = c2a@_FillValue
pha1(0,:) = pha1@_FillValue
;pha1 = where(c2a.lt.0.05, pha1@_FillValue, pha1) ; mask
n = 15
pha2 = STC(n,:,{-nWavePlt:nWavePlt})
pha2_at_long_name = "asymmetric phase-2"
pha2@_FillValue = c2a@_FillValue
pha2(0,:) = pha2@_FillValue
pha2 = where(c2a.lt.0.05, pha2@_FillValue, pha2) ; mask
if (opt .and. isatt(opt,"pltProb") ) then
np = ind(c2a_at_prob_coh2 .eq. opt_at_pltProb)
if (.not.ismissing(np)) then
c2a = where(c2a.lt.STC_at_prob_coh2(np), c2a@_FillValue, c2s)
pha1 = where(ismissing(c2a), pha1@_FillValue, pha1)
pha2 = where(ismissing(c2a), pha2@_FillValue, pha2)
end if
end if
if (opt .and. isatt(opt,"coh2Cutoff") ) then
c2a = where(c2a.lt.opt_at_coh2Cutoff, c2s@_FillValue, c2s)
end if
if (opt .and. isatt(opt,"phaseCutoff") ) then
pha1 = where(c2a.lt.opt_at_phaseCutoff, pha1@_FillValue, pha1)
pha2 = where(c2a.lt.opt_at_phaseCutoff, pha2@_FillValue, pha2)
end if
if (opt .and. isatt(opt,"elimIsoVals") .and. .not.opt_at_elimIsoVals) then
mjo_elimIsolatedValues (c2a, pha1, pha2, 2)
mjo_elimIsolatedValues (c2a, pha1, pha2, 1)
end if
if (plotPhase) then
scl_one = sqrt(1./(pha1^2 + pha2^2))
pha1 = scl_one*pha1
pha2 = scl_one*pha2
plot(1) = gsn_csm_vector_scalar(wks, pha1, pha2, c2a, res)
else
plot(1) = gsn_csm_contour(wks, c2a, res)
end if
duma = new (25, "graphic") ; for _add_ *a*symmetric graphical objects
plot(1) = addHorVertLinesCross(wks, plot(1), nWavePlt, duma)
dumdca = new ( 25, "graphic") ; dispersion curves asymmetric (dcs)
dumdca( 0) = gsn_add_polyline(wks,plot(1),Apzwn(0,0,:),Afreq(0,0,:),dcres)
dumdca( 1) = gsn_add_polyline(wks,plot(1),Apzwn(0,1,:),Afreq(0,1,:),dcres)
dumdca( 2) = gsn_add_polyline(wks,plot(1),Apzwn(0,2,:),Afreq(0,2,:),dcres)
dumdca( 3) = gsn_add_polyline(wks,plot(1),Apzwn(1,0,:),Afreq(1,0,:),dcres)
dumdca( 4) = gsn_add_polyline(wks,plot(1),Apzwn(1,1,:),Afreq(1,1,:),dcres)
dumdca( 5) = gsn_add_polyline(wks,plot(1),Apzwn(1,2,:),Afreq(1,2,:),dcres)
dumdca( 6) = gsn_add_polyline(wks,plot(1),Apzwn(2,0,:),Afreq(2,0,:),dcres)
dumdca( 7) = gsn_add_polyline(wks,plot(1),Apzwn(2,1,:),Afreq(2,1,:),dcres)
dumdca( 8) = gsn_add_polyline(wks,plot(1),Apzwn(2,2,:),Afreq(2,2,:),dcres)
dumdca(10) = gsn_add_text(wks,plot(1),"MRG",-10.0,.15,txres)
dumdca(11) = gsn_add_text(wks,plot(1),"n=2 IG",-3.0,.58,txres)
dumdca(12) = gsn_add_text(wks,plot(1),"n=0 EIG",6.5,.40,txres)
dumdca(13) = gsn_add_text(wks,plot(1),"h=50",-10.0,.78,txres)
dumdca(14) = gsn_add_text(wks,plot(1),"h=25",-10.0,.63,txres)
dumdca(15) = gsn_add_text(wks,plot(1),"h=12",-10.0,.51,txres)
resP = True
resP_at_gsnMaximize = True
resP_at_gsnPanelLabelBar = True
;resP_at_lbLabelAutoStride = True
resP_at_lbLabelStride = 2 ; every other one
;resP_at_lbLabelFontHeightF = 0.007 ; make labels smaller
resP_at_cnLabelBarEndLabelsOn= True
resP_at_cnLabelBarEndStyle = "ExcludeOuterBoxes"
resP_at_cnLevelSelectionMode = res_at_cnLevelSelectionMode
resP_at_cnMinLevelValF = res_at_cnMinLevelValF
resP_at_cnMaxLevelValF = res_at_cnMaxLevelValF ; correlation^2 = 0.8
resP_at_cnLevelSpacingF = res_at_cnLevelSpacingF
if (opt .and. isatt(opt,"txString") ) then
resP_at_txString = opt_at_txString
end if
if (opt .and. isatt(opt,"gsnPaperOrientation") ) then
resP_at_gsnPaperOrientation = opt_at_gsnPaperOrientation
end if
if (opt .and. isatt(opt,"lbOrientation") ) then
resP_at_lbOrientation = "vertical"
gsn_panel(wks,plot,(/2,1/),resP)
else
gsn_panel(wks,plot,(/1,2/),resP)
end if
if (pltType.eq."png") then
if (isatt(opt,"pltConvert")) then
pltConvert = opt_at_pltConvert ; convert options
else
pltConvert = "-rotate -90" ; default
end if
system("convert "+pltConvert+" "+pltPath+".eps "+pltPath+".png")
system("/bin/rm -f "+pltPath+".eps")
end if
end
; ---------------------
undef ("mjo_phase_background")
function mjo_phase_background (wks:graphic , opt[1]:logical )
begin
if (opt .and. isatt(opt,"axisExtent") .and. opt_at_axisExtent.gt.0) then
axisExtent = opt_at_axisExtent
else
axisExtent = 4
end if
nPhase = 8
res = True
res_at_gsnDraw = False
res_at_gsnFrame = False
;res_at_gsnPaperOrientation = "portrait"
res_at_vpXF = 0.1 ; default=0.2
res_at_vpYF = 0.8 ; default=0.8
res_at_trYMinF = -axisExtent
res_at_trYMaxF = axisExtent
res_at_trXMinF = -axisExtent
res_at_trXMaxF = axisExtent + 0.05
vpExtent = 0.45
res_at_vpWidthF = vpExtent ; default=0.6
res_at_vpHeightF = vpExtent ; default=0.6
res_at_tmXBFormat = "f" ; integer, no decimal
res_at_tmYLFormat = "f" ; 3 rather than 3.0
res_at_tmXBLabelDeltaF = -0.75 ; move label closer to tm
res_at_tmYLLabelDeltaF = -0.75
labelFontHeightF = 0.0167 ; default=0.0167
res_at_tiXAxisFontHeightF= labelFontHeightF
res_at_tiYAxisFontHeightF= labelFontHeightF
res_at_tiDeltaF = 1.25 ; default=1.5
res_at_xyLineThicknessF = 1
rad = 4.*atan(1.0)/180.
if (opt .and. isatt(opt,"radius") .and. opt_at_radius.gt.0) then
radius = opt_at_radius
else
radius= 1.0 ; unit length (default)
end if
nCirc = 361
xCirc = radius*cos(fspan(0, 360, nCirc)*rad)
yCirc = radius*sin(fspan(0, 360, nCirc)*rad)
xCirc_at_long_name = "Phase 2 (Indian Ocean) Phase 3"
yCirc_at_long_name = "Phase 1 (Western Hem, Africa) Phase 8"
;res_at_gsnCenterStringOrthogonalPosF =
res_at_gsnCenterStringFontHeightF = labelFontHeightF
res_at_gsnCenterString = "Phase 7 (Western Pacific) Phase 6"
if (opt .and. isatt(opt,"tiMainString")) then
res_at_tiMainString = opt_at_tiMainString
end if
plt = gsn_csm_xy (wks, xCirc, yCirc, res) ; base plot
txres = True
txres_at_txFontHeightF = labelFontHeightF
txres_at_txAngleF = -90
txid = gsn_create_text(wks, "Phase 5 (Maritime) Phase 4", txres)
amres = True
amres_at_amParallelPosF = 0.575
amres_at_amOrthogonalPosF = 0.0
amres_at_amJust = "CenterCenter"
ann3 = gsn_add_annotation(plt, txid, amres)
;************************************************
; Add the phase separator lines.
; Each line must be associated with a unique dummy variable.
; Draw each line separately. Each line must contain two
; points (begin/end).
;************************************************
plres = True ; polyline mods desired
;plres_at_gsLineColor = "red" ; color of lines
plres_at_gsLineThicknessF = 1.0 ; thickness of lines
if (opt .and. isatt(opt,"gsLineDashPattern")) then
plres_at_gsLineDashPattern = opt_at_gsLineDashPattern
else
plres_at_gsLineDashPattern = 1
end if
c45 = radius*cos(45*rad)
E = axisExtent ; convenience
R = radius
phaLine = new ( (/nPhase,4/), "float", "No_FillValue")
phaLine(0,:) = (/ R, E, 0, 0/) ; (xBegin,xEnd,yBegin,yEnd)
phaLine(1,:) = (/ c45, E, c45, E/)
phaLine(2,:) = (/ 0, 0, R, E/)
phaLine(3,:) = (/-c45, -E, c45, E/)
phaLine(4,:) = (/ -R, -E, 0, 0/)
phaLine(5,:) = (/-c45, -E,-c45, -E/)
phaLine(6,:) = (/ 0, 0, -R, -E/)
phaLine(7,:) = (/ c45, E,-c45, -E/)
do i =0,nPhase-1
plt@$unique_string("dum")$ = gsn_add_polyline(wks,plt \
,(/phaLine(i,0),phaLine(i,1)/) \
,(/phaLine(i,2),phaLine(i,3)/) ,plres)
end do
return(plt)
end
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Tue Apr 21 2009 - 09:44:55 MDT
This archive was generated by hypermail 2.2.0 : Tue Apr 21 2009 - 11:04:41 MDT