[ncl-talk] fatal:divide: Division by 0, Can't continuefatal:Div: operator failed, can't continue fatal:["Execute.c":8635]:Execute: Error occurred at or near line 129
Dennis Shea
shea at ucar.edu
Mon Dec 23 08:32:22 MST 2019
No!
---
mjoclivar_14 *requires multiple years of daily anomaly data*. Note the
example uses anomalies from:
ymdStrt = 19800101 ; start yyyymmdd
ymdLast = 20051231 ; last
Note the Lanczos band-pass filter specified by the MJO Clivar
diagnostics. See top of section.
;************************************************
; create BandPass Filter
;************************************************
ihp = 2 ; bpf=>band pass filter
*nWgt = 201*
sigma = 1.0 ; Lanczos sigma
fca = 1./100.
fcb = 1./20.
wgt = *filwgts_lanczos*
<http://www.ncl.ucar.edu/Document/Functions/Built-in/filwgts_lanczos.shtml>
(nWgt, ihp, fca, fcb, sigma )
*Hence to capture year 2018 you need (nWgt/2=201/2=) 150 or more days
into year 2019.*
---
*After* doing mjoclivar_14 for multiple years, you then specifya
period in mjoclivar_15.
On Fri, Dec 20, 2019 at 7:23 PM Rahpeni Fajarianti <
rahpenifajarianti at gmail.com> wrote:
> Thank you.
>
> Now, I use long term means for daily values for 1981-2010 as base
> climatology data. It works.
> I got min max values for :
> - OLR anomaly (min=-185.5106410148043 max=77.65257890949087)
> - 200hPa zonal wind anomaly (min=-48.11269687686956
> max=57.21568789939967)
> - 850hPa zonal wind anomaly (min=-31.22578327383305
> max=41.24650545117836).
>
> then, I use mjoclivar_14.ncl to get multivariate eof for data 2018. I
> attach my script.
>
> Some error shown like this:
> sub drveof: ier,jopt= 10 0 returned from vcmssm/crmssm
> warning:eofunc: 10 eigenvectors failed to converge
> (0) ==============
> Variable: eof_cdata
> Type: double
> Total Size: 6960 bytes
> 870 values
> Number of Dimensions: 2
> Dimensions and sizes: [2] x [435]
> Coordinates:
> Number Of Attributes: 5
> _FillValue : -999000000
> method : no transpose
> matrix : covariance
> pcvar : ( -9.99e+08, -9.99e+08 )
> eval : ( -999000000, -999000000 )
> (0)
> (0) min=-999000000 max=-999000000
> warning:eofunc_ts: 2 eigenvectors failed to converge
> (0) ==============
>
> Variable: eof_ts_cdata
> Type: double
> Total Size: 5840 bytes
> 730 values
> Number of Dimensions: 2
> Dimensions and sizes: [2] x [365]
> Coordinates:
> Number Of Attributes: 3
> _FillValue : -999000000
> matrix : covariance
> ts_mean : ( -999000000, -999000000 )
> (0)
> (0) min=-999000000 max=-999000000
> warning:eofunc_ts: 2 eigenvectors failed to converge
> warning:eofunc_ts: 2 eigenvectors failed to converge
> warning:eofunc_ts: 2 eigenvectors failed to converge
> warning:avg: Entire input array contains missing values, can't compute
> average
>
> fatal:Subscript out of range, error in subscript #0
> fatal:An error occurred reading unknown
> fatal:["Execute.c":8635]:Execute: Error occurred at or near line 207
>
> this is the error line :
> ncl 207> lonmax_eof1 = ceof&lon(imax_olr_eof1) ; longitude of max
> value (i.e. suppressed convection)
> ncl 208> lonmax_eof2 = ceof&lon(imax_olr_eof2)
>
> The data I used:
> olranom.nc
> <https://drive.google.com/file/d/1xbmRiRO9VrTQ6qqKTfTq_rraiRSZFOXg/view?usp=drive_web>
> uanom.200.nc
> <https://drive.google.com/file/d/1u1IYOGmjLuWAR4izwLkAgEMzWDxrkZlX/view?usp=drive_web>
> uanom.850.nc
> <https://drive.google.com/file/d/1ZTtEukGralmUuwBWLGJ9ei6-N0yf199b/view?usp=drive_web>
> The script I used:
>
> On Sat, Dec 21, 2019 at 2:21 AM Dennis Shea <shea at ucar.edu> wrote:
>
>> [1]
>> PLEASE do not include your line-by-line script in you emails. It makes
>> the emails too long. Attach the script.
>>
>> [2]
>> As noted in mjoclivar_2.ncl, you must use *multiple years* of data to
>> create a base climatology. Specifically: "This example uses 26-years
>> (1980-2005) of daily data". ANomalies are derived as difference from some
>> base (multi-year) climatology.
>>
>> [3]
>> I speculate you are using daily values from 2018 only. That is ***NOT***
>> a climatology. [Look up definition of climatology.]
>> There are no anomalies when you use one year of data That is why
>> min-anomaly=0.0 md max-anomaly=0.0
>>
>>
>>
>>
>>
>>
>> On Fri, Dec 20, 2019 at 10:14 AM Rahpeni Fajarianti via ncl-talk <
>> ncl-talk at ucar.edu> wrote:
>>
>>> I did it. Thanks.
>>>
>>> I have been trying many times with this data to get olr and zonal wind
>>> anomaly and finally I got both min and max data values.
>>>
>>> *This is the result when I use printMinMax statement :*
>>>
>>> Variable: olr
>>> Type: double
>>> Total Size: 5504200 bytes
>>> 688025 values
>>> Number of Dimensions: 3
>>> Dimensions and sizes: [time | 365] x [lat | 13] x [lon | 145]
>>> Coordinates:
>>> time: [ 0..524160]
>>> lat: [ -15.. 15]
>>> lon: [ 0.. 360]
>>> Number Of Attributes: 1
>>> _FillValue : -999000000
>>> ncl 19> printMinMax(olr,0)
>>> *(0) min=-185.5106410148043 max=77.65257890949087*
>>> ncl 20> print("------------------------------")
>>> (0) ------------------------------
>>> ncl 21> printVarSummary(u200)
>>>
>>> Variable: u200
>>> Type: double
>>> Total Size: 5504200 bytes
>>> 688025 values
>>> Number of Dimensions: 3
>>> Dimensions and sizes: [time | 365] x [lat | 13] x [lon | 145]
>>> Coordinates:
>>> time: [ 0..524160]
>>> lat: [ -15.. 15]
>>> lon: [ 0.. 360]
>>> Number Of Attributes: 1
>>> _FillValue : -999000000
>>> ncl 22> printMinMax(u200,0)
>>> *(0) min=-48.11269687686956 max=57.21568789939967*
>>> ncl 23> print("------------------------------")
>>> (0) ------------------------------
>>> ncl 24> printVarSummary(u850)
>>>
>>> Variable: u850
>>> Type: double
>>> Total Size: 5504200 bytes
>>> 688025 values
>>> Number of Dimensions: 3
>>> Dimensions and sizes: [time | 365] x [lat | 13] x [lon | 145]
>>> Coordinates:
>>> time: [ 0..524160]
>>> lat: [ -15.. 15]
>>> lon: [ 0.. 360]
>>> Number Of Attributes: 1
>>> _FillValue : -999000000
>>> ncl 25> printMinMax(u850,0)
>>> *(0) min=-31.22578327383305 max=41.24650545117836*
>>>
>>> I want to run mjoclivar_14.ncl (
>>> https://www.ncl.ucar.edu/Applications/mjoclivar.shtml), but I got is
>>> always both min and max data value -999000000. Could somebody help me
>>> making multivariate EOF with this file :
>>>
>>> OLR anomaly data:
>>> olranom.nc
>>> <https://drive.google.com/file/d/1bH67Flh4XgfYfSfsrBRPd0F7Gwq9qce5/view?usp=drive_web>
>>> 200 hPa zonal wind anomaly data :
>>> uanom.200.nc
>>> <https://drive.google.com/file/d/1Y_JmVgDatlxeQ2VsPSD4g-dh6Cc0XNna/view?usp=drive_web>
>>> 850 hPa zonal wind anomaly data:
>>> uanom.850.nc
>>> <https://drive.google.com/file/d/1TwaQL7zF8J3wIyCs8mmnQC1lqimPFsCV/view?usp=drive_web>
>>>
>>> This is my script for mjoclivar_14 :
>>> ncl 0> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/diagnostics_cam.ncl"
>>> dim_rmvmean_n( olr, 0) ; (time,lon)
>>> u850 = dimncl 1>
>>> ;******************************************************
>>> ncl 2> ;
>>> ncl 3> ; mjoclivar_14.ncl
>>> ncl 4> ;
>>> ncl 5> ;***********************************************************
>>> ncl 6> ; Combined EOFs
>>> ncl 7> ; Latest Update: July, 2016: Eun-Pa Lim; Bureau of Meteorology,
>>> Australia
>>> ncl 8> ;***********************************************************
>>> ncl 9> ;;
>>> ncl 10> ;; The following are automatically loaded from 6.2.0 onward
>>> ncl 11> ;;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
>>> ncl 12> ;;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
>>> ncl 13> ;;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
>>> ncl 14>
>>> ncl 15> undef("read_rename")
>>> ncl 16> function read_rename(f[1]:file, varName[1]:string \
>>> var_olr = dim_avg_n_Wrap( var_olr , 0)
>>> zavg_var_u850 = dim_avg_n_Wrap( var_u850, 0)
>>> zavg_var_u200 = dim_avg_n_Wrap(ncl 16>
>>> ,iStrt[1]:integer, iLast[1]:integer \
>>> ncl 16> ,latS[1]:numeric , latN[1]:numeric )
>>> ncl 17> ; Utility to force specific named dimensions
>>> ncl 18> ; This is done for historical reasons (convenience)
>>> ncl 19> begin
>>> ncl 20> work = f->$varName$(iStrt:iLast,{latS:latN},:) ;
>>> (time,lat,lon)
>>> ncl 21> work!0 = "time" ; CAM
>>> model names
>>> ncl 22> work!1 = "lat"
>>> ncl 23> work!2 = "lon"
>>> ncl 24> return(work)
>>> ncl 25> end
>>> ncl 26> ; =========================> MAIN
>>> <==============================
>>> ncl 27> begin
>>> ncl 28> neof = 2
>>> ncl 29>
>>> ncl 30> latS = -15
>>> ncl 31> latN = 15
>>> ncl 32>
>>> mncl 33> ymdStrt = 20180101 ; start yyyymmdd
>>> ncl 34> ymdLast = 20181231 ; last
>>> ncl 35>
>>> yncl 36> yrStrt = ymdStrt/10000
>>> ncl 37> yrLast = ymdLast/10000
>>> ncl 38>
>>> ancl 39> pltDir = "./" ; plot directory
>>> ncl 40> pltType = "png"
>>> ncl 41> pltName = "mjoclivar"
>>> ncl 42>
>>> ncl 43> diri = "./" ; input directory
>>> ncl 44>
>>> oncl 45> filolr = "olranom.nc"
>>> ncl 46> filu200 = "uanom.200.nc"
>>> ncl 47> filu850 = "uanom.850.nc"
>>> ncl 48>
>>> ncl 49> ;************************************************
>>> ncl 50> ; create BandPass Filter
>>> ncl 51> ;************************************************
>>> oncl 52> ihp = 2 ; bpf=>band pass
>>> filter
>>> ncl 53> nWgt = 201
>>> ncl 54> sigma = 1.0 ; Lanczos sigma
>>> ncl 55> fca = 1./100.
>>> ncl 56> fcb = 1./20.
>>> ncl 57> wgt = filwgts_lanczos (nWgt, ihp, fca, fcb, sigma )
>>> ncl 58>
>>> *ncl 59> ;***********************************************************
>>> ncl 60> ; Find the indices corresponding to the start/end times
>>> ncl 61> ;***********************************************************
>>> ncl 62> f = addfile (diri+filolr , "r")
>>> ncl 63> TIME = f->time ; days since ...
>>> ncl 64>
>>> ncl 65> YMD = cd_calendar(TIME, -2) ; entire (time,6)
>>> ncl 66>
>>> ncl 67> iStrt = ind(YMD.eq.ymdStrt) ; index start
>>> ncl 68> iLast = ind(YMD.eq.ymdLast) ; index last
>>> ncl 69> delete([/ TIME, YMD /])
>>> ncl 70>
>>> ncl 71> ;***********************************************************
>>> ncl 72> ; Read anomalies
>>> ncl 73> ;***********************************************************
>>> ncl 74>
>>> encl 75> work = read_rename(f,"olranom",iStrt,iLast,latS,latN) ;
>>> (time,lat,lon)
>>> ncl 76> OLR = dim_avg_n_Wrap(work, 1) ;
>>> (time,lon)
>>> ncl 77> delete(work)
>>> ncl 78>
>>> ncl 79> f = addfile (diri+filu850 , "r")
>>>
>>> ncl 80> work = read_rename(f,"uanom",iStrt,iLast,latS,latN) ;
>>> (time,lat,lon)
>>> ncl 81> U850 = dim_avg_n_Wrap(work, 1) ; (time,lon)
>>> ncl 82> delete(work)
>>> ncl 83>
>>> ncl 84> f = addfile (diri+filu200 , "r")
>>>
>>> ncl 85> work = read_rename(f,"uanom",iStrt,iLast,latS,latN) ;
>>> (time,lat,lon)
>>> ncl 86> U200 = dim_avg_n_Wrap(work, 1) ; (time,lon)
>>> ncl 87>
>>> 0ncl 88> dimw = dimsizes( work )
>>> ncl 89> ntim = dimw(0)
>>> ncl 90> nlat = dimw(1)
>>> )ncl 91> mlon = dimw(2)
>>> uncl 92> delete(work)
>>> ncl 93>
>>> oncl 94> lon = OLR&lon
>>> ncl 95> time = OLR&time
>>> ncl 96> date = cd_calendar(time, -2) ; yyyymmdd
>>> ncl 97>
>>> encl 98> ;************************************************
>>> ncl 99> ; Apply the band pass filter to the original anomalies
>>> ncl 100> ;************************************************
>>> ncl 101> olr = wgt_runave_n_Wrap ( OLR, wgt, 0, 0) ; (time,lon)
>>> ncl 102> u850 = wgt_runave_n_Wrap (U850, wgt, 0, 0)
>>> ncl 103> u200 = wgt_runave_n_Wrap (U200, wgt, 0, 0)
>>> ncl 104>
>>> *ncl 105> ;************************************************
>>> ncl 106> ; remove temporal means of band pass series: *not* necessary
>>> ncl 107> ;************************************************
>>> ncl 108> olr = dim_rmvmean_n( olr, 0) ; (time,lon)
>>> ncl 109> u850 = dim_rmvmean_n(u850, 0)
>>> ncl 110> u200 = dim_rmvmean_n(u200, 0)
>>> ncl 111>
>>> nncl 112> ;************************************************
>>> ncl 113> ; Compute the temporal variance at each lon
>>> ncl 114> ;************************************************
>>> ncl 115> var_olr = dim_variance_n_Wrap( olr, 0) ; (lon)
>>> ncl 116> var_u850 = dim_variance_n_Wrap(u850, 0)
>>> ncl 117> var_u200 = dim_variance_n_Wrap(u200, 0)
>>> ncl 118>
>>> ncl 119> ;************************************************
>>> ncl 120> ; Compute the zonal mean of the temporal variance
>>> ncl 121> ;************************************************
>>> ncl 122> zavg_var_olr = dim_avg_n_Wrap( var_olr , 0)
>>> ncl 123> zavg_var_u850 = dim_avg_n_Wrap( var_u850, 0)
>>> ncl 124> zavg_var_u200 = dim_avg_n_Wrap( var_u200, 0)
>>> ncl 125>
>>> rncl 126> ;************************************************
>>> ancl 127> ; Normalize by sqrt(avg_var*)
>>> ncl 128> ;************************************************
>>> ncl 129> zsqrt =
>>> where(zavg_var_olr.gt.0,sqrt(zavg_var_olr),zavg_var_olr at _FillValue)
>>> ncl 130> olr = olr/zsqrt ; (time,lon)
>>> ;ncl 131> z850sqrt
>>> =where(zavg_var_u850.gt.0,sqrt(zavg_var_u850),zavg_var_u850 at _FillValue)
>>> ncl 132> u850 = u850/z850sqrt
>>> ncl 133> z200sqrt=
>>> where(zavg_var_u200.gt.0,sqrt(zavg_var_u200),zavg_var_u200 at _FillValue)
>>> )ncl 134> u200 = u200/z200sqrt
>>> ncl 135>
>>> cncl 136> ;************************************************
>>> ncl 137> ; Combine the normalized data into one variable
>>> ncl 138> ;************************************************
>>> ncl 139> cdata = new ( (/3*mlon,ntim/), typeof(olr),
>>> getFillValue(olr))
>>> ncl 140> do ml=0,mlon-1
>>> ncl 141> cdata(ml ,:) = (/ olr(:,ml) /)
>>> ncl 142> cdata(ml+ mlon,:) = (/ u850(:,ml) /)
>>> ncl 143> cdata(ml+2*mlon,:) = (/ u200(:,ml) /)
>>> ncl 144> end do
>>> ncl 145>
>>> =ncl 146> ;************************************************
>>> ncl 147> ; Compute **combined** EOF; Sign of EOF is arbitrary
>>> ncl 148> ;************************************************
>>> ncl 149> eof_cdata = eofunc(cdata , neof, False) ;
>>> (neof,3*mlon)
>>> ncl 150> print("==============")
>>> ncl 151> printVarSummary(eof_cdata)
>>> ncl 152> printMinMax(eof_cdata, True)
>>> ncl 153>
>>> )ncl 154> eof_ts_cdata = eofunc_ts(cdata,eof_cdata,False) ;
>>> (neof,3*ntim)
>>> ncl 155> print("==============")
>>> ncl 156> printVarSummary(eof_ts_cdata)
>>> ncl 157> printMinMax(eof_ts_cdata, True)
>>> ncl 158>
>>> ncl 159> ;************************************************
>>> ncl 160> ; For clarity, explicitly extract each variable. Create time
>>> series
>>> ncl 161> ;************************************************
>>> ncl 162>
>>> tncl 163> nvar = 3 ; "olr", "u850", "u200"
>>> ncl 164> ceof = new( (/nvar,neof,mlon/), typeof(cdata),
>>> getFillValue(cdata))
>>> ncl 165>
>>> dncl 166> do n=0,neof-1
>>> ncl 167> ceof(0,n,:) = eof_cdata(n,0:mlon-1) ; olr
>>> ncl 168> ceof(1,n,:) = eof_cdata(n,mlon:2*mlon-1) ; u850
>>> ncl 169> ceof(2,n,:) = eof_cdata(n,2*mlon:) ; u200
>>> ncl 170> end do
>>> ncl 171>
>>> (ncl 172> ceof!0 = "var"
>>> ncl 173> ceof!1 = "eof"
>>> ncl 174> ceof!2 = "lon"
>>> ncl 175> ceof&lon = olr&lon
>>> ncl 176>
>>> incl 177> ceof_ts = new( (/nvar,neof,ntim/), typeof(cdata),
>>> getFillValue(cdata))
>>> ncl 178> ceof_ts(0,:,:) = eofunc_ts_Wrap(
>>> olr(lon|:,time|:),ceof(0,:,:),False) ; (0,neof,ntim)
>>> ncl 179> ceof_ts(1,:,:) =
>>> eofunc_ts_Wrap(u850(lon|:,time|:),ceof(1,:,:),False) ; (1,neof,ntim)
>>> ncl 180> ceof_ts(2,:,:) =
>>> eofunc_ts_Wrap(u200(lon|:,time|:),ceof(2,:,:),False) ; (2,neof,ntim)
>>> ncl 181>
>>> tncl 182> ;**********************************************t*
>>> ncl 183> ; Add code contributed by Marcus N. Morgan, Florida Institute
>>> of Technology; Feb 2015
>>> ncl 184> ; Calculate % variance (pcv_ )accounted for by OLR, U850 and
>>> U200
>>> ncl 185> ;************************************************
>>> ncl 186>
>>> ncl 187> pcv_eof_olr = new(neof,typeof(ceof))
>>> ncl 188> pcv_eof_u850 = new(neof,typeof(ceof))
>>> ncl 189> pcv_eof_u200 = new(neof,typeof(ceof))
>>> ncl 190>
>>> ncl 191> do n=0,neof-1
>>> ncl 192> pcv_eof_olr(n) = avg((ceof(0,n,:)*sqrt(ceof at eval
>>> (n)))^2)*100
>>> ncl 193> pcv_eof_u850(n) = avg((ceof(1,n,:)*sqrt(ceof at eval
>>> (n)))^2)*100
>>> ncl 194> pcv_eof_u200(n) = avg((ceof(2,n,:)*sqrt(ceof at eval
>>> (n)))^2)*100
>>> ncl 195> ;;print("pcv: neof="+(n+1)+": "+pcv_eof_olr(n)+"
>>> "+pcv_eof_u850(n)+" "+pcv_eof_u200(n))
>>> ncl 196> end do
>>> ncl 197>
>>> ncl 198> ;************************************************
>>> ncl 199> ; Change sign of EOFs for spatial structures of PC1 and PC2
>>> ncl 200> ; to represent convection over the tropical Indian Ocean and
>>> the tropical western Pacific Ocean, respectively
>>> ncl 201> ; (Ad hoc approach)
>>> ncl 202> ;************************************************
>>> ncl 203>
>>> ncl 204> imax_olr_eof1 = maxind(ceof(0,0,:))
>>> ncl 205> imax_olr_eof2 = maxind(ceof(0,1,:))
>>> ncl 206>
>>> )ncl 207> lonmax_eof1 = ceof&lon(imax_olr_eof1) ; longitude of
>>> max value (i.e. suppressed convection)
>>> ncl 208> lonmax_eof2 = ceof&lon(imax_olr_eof2)
>>> ncl 209>
>>> sncl 210> if (lonmax_eof1.ge.60 .and. lonmax_eof1.lt.180) then ;
>>> Change the sign of EOF1
>>> ncl 211> ceof(:,0,:) = -ceof(:,0,:) ; if
>>> OLR is positive
>>> ncl 212> ceof_ts(:,0,:) = -ceof_ts(:,0,:) ; over
>>> the tropical Indian Ocean
>>> ncl 213> eof_cdata(0,:) = -eof_cdata(0,:)
>>> ncl 214> eof_ts_cdata(0,:) = -eof_ts_cdata(0,:)
>>> ncl 215> end if
>>> ncl 216>
>>> encl 217> if (lonmax_eof2.ge.120 .and. lonmax_eof2.lt.180) then ;
>>> Change the sign of EOF2
>>> ncl 218> ceof(:,1,:) = -ceof(:,1,:) ; if
>>> OLR is positive
>>> ncl 219> ceof_ts(:,1,:) = -ceof_ts(:,1,:) ; over
>>> the tropical western Pacific Ocean
>>> ncl 220> eof_cdata(1,:) = -eof_cdata(1,:)
>>> ncl 221> eof_ts_cdata(1,:) = -eof_ts_cdata(1,:)
>>> ncl 222> end if
>>> ncl 223>
>>> ncl 224> print("==============")
>>> ncl 225> printVarSummary(eof_cdata)
>>> ncl 226> printMinMax(eof_cdata, True)
>>> nncl 227>
>>> ncl 228> ;************************************************
>>> ncl 229> ; Compute cross correlation of each variable's EOF time series
>>> at zero-lag
>>> ncl 230> ;************************************************
>>> ncl 231> r_olr_u850 = escorc(ceof_ts(0,:,:) , ceof_ts(1,:,:) ) ;
>>> (neof)
>>> ncl 232> r_olr_u200 = escorc(ceof_ts(0,:,:) , ceof_ts(2,:,:) )
>>> ncl 233> r_u850_u200 = escorc(ceof_ts(1,:,:) , ceof_ts(2,:,:) )
>>> ncl 234>
>>> ncl 235> print("==============")
>>> ncl 236> do n=0,neof-1
>>> ncl 237> print("neof="+n \
>>> ncl 237> +" r_olr_u850=" +sprintf("%4.3f",r_olr_u850(n)) \
>>> ncl 237> +" r_olr_u200=" +sprintf("%4.3f",r_olr_u200(n)) \
>>> ncl 237> +" r_u850_u200="+sprintf("%4.3f",r_u850_u200(n)) )
>>> ncl 238> end do
>>> ncl 239> print("==============")
>>> ncl 240>
>>> ncl 241> ;************************************************
>>> ncl 242> ; Compute cross correlation of the multivariate EOF; EOF 1 vs
>>> EOF 2
>>> ncl 243> ;************************************************
>>> ncl 244>
>>> ncl 245> mxlag = 25
>>> ncl 246> rlag_01 = esccr(eof_ts_cdata(0,:),eof_ts_cdata(1,:), mxlag)
>>> ; (N,mxlag+1)
>>> ncl 247> rlag_10 = esccr(eof_ts_cdata(1,:),eof_ts_cdata(0,:), mxlag)
>>> ; (N,mxlag+1)
>>> ncl 248> ccr_12 = new ( (/2*mxlag+1/), float)
>>> ncl 249>
>>> ncl 250> ccr_12(mxlag:) = rlag_10(0:mxlag)
>>> ncl 251> ccr_12(0:mxlag) = rlag_01(::-1) ; reverse order
>>> ncl 252> ;;print(ccr_12)
>>> ncl 253>
>>> ncl 254>
>>> ncl 255> ;************************************************
>>> ncl 256> ; Normalize the multivariate EOF 1&2 component time series
>>> ncl 257> ; Compute (PC1^2+PC2^2): values > 1 indicate "strong" periods
>>> encl 258> ;************************************************
>>> ncl 259> eof_ts_cdata(0,:) =
>>> eof_ts_cdata(0,:)/stddev(eof_ts_cdata(0,:))
>>> ncl 260> eof_ts_cdata(1,:) =
>>> eof_ts_cdata(1,:)/stddev(eof_ts_cdata(1,:))
>>> ncl 261>
>>> ncl 262> mjo_ts_index = eof_ts_cdata(0,:)^2 + eof_ts_cdata(1,:)^2
>>> ncl 263> mjo_ts_index_smt = runave(mjo_ts_index, 91, 0) ; 91-day
>>> running mean
>>> ncl 264>
>>> ncl 265> nGood = num(.not.ismissing(mjo_ts_index)) ; #
>>> non-missing
>>> ncl 266> nStrong = num(mjo_ts_index .ge. 1.0)
>>> ncl 267> print("nGood="+nGood+" nStrong="+nStrong+"
>>> nOther="+(nGood-nStrong))
>>> ncl 268>
>>> rncl 269> ;************************************************
>>> uncl 270> ; Write PC results to netCDF for use in another example.
>>> ncl 271> ;************************************************
>>> ncl 272> mjo_ts_index!0 = "time"
>>> t
>>> ncl 273> mjo_ts_index&time = time
>>> ncl 274> mjo_ts_index at long_name = "MJO PC INDEX"
>>>
>>>
>>> ;****ncl 275> mjo_ts_index at info = "(PC1^2 + PC2^2)"
>>> ncl 276>
>>> ncl 277> PC1 = eof_ts_cdata(0,:)
>>> ncl 278> PC1!0 = "time"
>>> ********************************************
>>> ; MJO "strong" index
>>> ;********************ncl 279> PC1&time = time
>>> ncl 280> PC1 at long_name = "PC1"
>>> ncl 281> PC1 at info = "PC1/stddev(PC1)"
>>> ncl 282>
>>> oncl 283> PC2 = eof_ts_cdata(1,:)
>>> ncl 284> PC2!0 = "time"
>>> ncl 285> PC2&time = time
>>> ncl 286> PC2 at long_name = "PC2"
>>> ncl 287> PC2 at info = "PC2/stddev(PC2)"
>>> ncl 288>
>>> ncl 289> diro = "./"
>>> ncl 290> filo = "MJO_PC_INDEX.nc"
>>> ncl 291> system("/bin/rm -f "+diro+filo) ; remove any pre-existing
>>> file
>>> ncl 292> ncdf = addfile(diro+filo,"c") ; open output netCDF file
>>> ncl 293> ; make time an UNLIMITED
>>> dimension
>>> ncl 294> filedimdef(ncdf,"time",-1,True) ; recommended for most
>>> applications
>>> ncl 295> ; output variables directly
>>> ncl 296> ncdf->MJO_INDEX = mjo_ts_index
>>> ncl 297> ncdf->PC1 = PC1
>>> ncl 298> ncdf->PC2 = PC2
>>> ncl 299>
>>> ncl 300> ;------------------------------------------------------------
>>> ncl 301> ; PLOTS
>>> ncl 302> ;------------------------------------------------------------
>>> Pncl 303>
>>> ncl 304> yyyymmdd = cd_calendar(time, -2)
>>> ncl 305> yrfrac = yyyymmdd_to_yyyyfrac(yyyymmdd, 0.0)
>>> ncl 306> delete([/ yrfrac at long_name, lon at long_name /])
>>> ncl 307>
>>> ncl 308> day = ispan(-mxlag, mxlag, 1)
>>> ncl 309> ;day at long_name = "lag (day)"
>>> ncl 310>
>>> ncl 311> pltPath = pltDir+pltName
>>> ncl 312>
>>> ncl 313> wks = gsn_open_wks(pltType,pltPath)
>>> ncl 314> gsn_define_colormap(wks,"default")
>>> ncl 315> plot = new(3,graphic)
>>> ncl 316>
>>> ncl 317> ;************************************************
>>> ncl 318> ; Multivariate EOF plots
>>> ncl 319> ;************************************************
>>> ncl 320> rts = True
>>> ncl 321> rts at gsnDraw = False ; don't draw yet
>>> ncl 322> rts at gsnFrame = False ; don't advance frame yet
>>> ncl 323> rts at gsnScale = True ; force text scaling
>>>
>>> ncl 324>
>>> ncl 325> rts at vpHeightF = 0.40 ; Changes the aspect ratio
>>> ncl 326> rts at vpWidthF = 0.85
>>> ncl 327> rts at vpXF = 0.10 ; change start locations
>>> ncl 328> rts at vpYF = 0.75 ; the plot
>>> ncl 329> rts at xyLineThicknesses = (/2, 2, 2/)
>>> ncl 330> rts at xyLineColors = (/"black","red","blue"/)
>>> ncl 331> rts at gsnYRefLine = 0. ; reference line
>>>
>>> ncl 332> rts at trXMaxF = max(lon)
>>> ncl 333> rts at trXMinF = min(lon)
>>> ncl 334>
>>> ncl 335> rts at pmLegendDisplayMode = "Always" ; turn on
>>> legend
>>> ncl 336> rts at pmLegendSide = "Top" ; Change
>>> location of
>>> ncl 337> rts at pmLegendParallelPosF = 1.16 ; move
>>> units right
>>> ncl 338> rts at pmLegendOrthogonalPosF = -0.50 ; move
>>> units down
>>> ncl 339> rts at pmLegendWidthF = 0.15 ; Change
>>> width and
>>> ncl 340> rts at pmLegendHeightF = 0.15 ; height of
>>> legend.
>>> ncl 341> rts at lgLabelFontHeightF = 0.0175
>>> ncl 342>
>>> ncl 343>
>>> ncl 344> rtsP = True ; modify the
>>> panel plot
>>> ncl 345> ; rtsP at gsnMaximize = True ; large
>>> format
>>> ncl 346> rtsP at gsnPanelMainString = "Multivariate EOF: 15S-15N:
>>> "+yrStrt+"-"+yrLast
>>> ncl 347>
>>> ncl 348> do n=0,neof-1
>>> ncl 349> rts at xyExplicitLegendLabels = (/"OLR: "+sprintf("%4.1f",
>>> pcv_eof_u200(n)) +"%" \
>>> ncl 349> ,"U850: "+sprintf("%4.1f",
>>> pcv_eof_u850(n))+"%" \
>>> ncl 349> ,"U200: "+sprintf("%4.1f",
>>> pcv_eof_olr(n))+"%" /)
>>> ncl 350> rts at gsnLeftString = "EOF "+(n+1)
>>> ncl 351> rts at gsnRightString = sprintf("%3.1f",ceof at pcvar(n)) +"%"
>>> ncl 352> plot(n) = gsn_csm_xy (wks,lon,ceof(:,n,:),rts)
>>> ncl 353> end do
>>> ncl 354> gsn_panel(wks,plot(0:1),(/2,1/),rtsP) ; now draw as one
>>> plot
>>> ncl 355>
>>> ncl 356> ;-----------------------------------------
>>> ncl 357> ; The following doesn't work with some older versions of NCL
>>> ncl 358> ; With old versions, the user must delete each individually.
>>> ncl 359> ;-----------------------------------------
>>> ncl 360> delete([/ rts at xyExplicitLegendLabels, rts at pmLegendDisplayMode
>>> \
>>> ncl 360> , rts at xyLineThicknesses , rts at gsnLeftString
>>> \
>>> ncl 360> , rts at gsnRightString , rts at xyLineColors
>>> \
>>> ncl 360> , rts at trXMaxF , rts at trXMinF
>>> /] )
>>> ncl 361>
>>> ncl 362> lag = ispan(-mxlag,mxlag,1)
>>> ncl 363> lag at long_name = "lag (days)"
>>> ncl 364>
>>> ncl 365> plot(0) = gsn_csm_xy (wks, lag ,ccr_12,rts)
>>> ncl 366> rtsP at gsnPanelMainString = "Cross Correlation:
>>> Multivariate EOF: 15S-15N: " \
>>> ncl 366> + yrStrt+"-"+yrLast
>>> ncl 367> rtsP at gsnPaperOrientation = "portrait" ; force
>>> portrait
>>> ncl 368> gsn_panel(wks,plot(0),(/1,1/),rtsP) ; now draw as one plot
>>> ncl 369>
>>> ncl 370> ;************************************************
>>> ncl 371> ; MJO "strong" index
>>> ncl 372> ;************************************************
>>> ncl 373> rts at gsnYRefLine = 1.0
>>> ncl 374> rts at gsnYRefLineColor = "black"
>>> ncl 375> rts at xyMonoDashPattern = True
>>> ncl 376> rts at xyLineColors = (/"black", "blue"/)
>>> ncl 377> rts at xyLineThicknesses = (/1, 2/)
>>> ncl 378> rts at pmLegendDisplayMode = "Always" ; turn on
>>> legend
>>> ncl 379> rts at pmLegendWidthF = 0.12 ; Change
>>> width and
>>> ncl 380> rts at pmLegendHeightF = 0.10 ; height of
>>> legend.
>>> ncl 381> rts at pmLegendParallelPosF = 0.86 ; move
>>> units right
>>> ncl 382> rts at pmLegendOrthogonalPosF = -0.40 ; move
>>> units down
>>> ncl 383> rts at xyExplicitLegendLabels = (/"daily", "91-day runavg" /)
>>> ncl 384>
>>> ncl 385> mjo_ind_plt = new ( (/2,ntim/), typeof(mjo_ts_index))
>>> ncl 386> mjo_ind_plt(0,:) = mjo_ts_index
>>> ncl 387> mjo_ind_plt(1,:) = (/ mjo_ts_index_smt /)
>>> ncl 388> plot(0) = gsn_csm_xy(wks, yrfrac,mjo_ind_plt,rts)
>>> ncl 389>
>>> ncl 390> rtsP at gsnPanelMainString = "MJO Index: (PC1^2+ PC2^2) :
>>> 15S-15N: "+yrStrt+"-"+yrLast
>>> ncl 391> gsn_panel(wks,plot(0),(/1,1/),rtsP) ; now draw as one plot
>>> ncl 392>
>>> ncl 393> end
>>>
>>>
>>> sub drveof: ier,jopt= 10 0 returned from vcmssm/crmssm
>>> warning:eofunc: 10 eigenvectors failed to converge
>>> (0) ==============
>>>
>>> Variable: eof_cdata
>>> Type: double
>>> Total Size: 6960 bytes
>>> 870 values
>>> Number of Dimensions: 2
>>> Dimensions and sizes: [2] x [435]
>>> Coordinates:
>>> Number Of Attributes: 5
>>> _FillValue : -999000000
>>> method : no transpose
>>> matrix : covariance
>>> pcvar : ( -9.99e+08, -9.99e+08 )
>>> eval : ( -999000000, -999000000 )
>>> (0)
>>> (0) min=-999000000 max=-999000000
>>> warning:eofunc_ts: 2 eigenvectors failed to converge
>>> (0) ==============
>>>
>>> Variable: eof_ts_cdata
>>> Type: double
>>> Total Size: 5840 bytes
>>> 730 values
>>> Number of Dimensions: 2
>>> Dimensions and sizes: [2] x [365]
>>> Coordinates:
>>> Number Of Attributes: 3
>>> _FillValue : -999000000
>>> matrix : covariance
>>> ts_mean : ( -999000000, -999000000 )
>>> (0)
>>> (0) min=-999000000 max=-999000000
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>> *warning:eofunc_ts: 2 eigenvectors failed to convergewarning:eofunc_ts:
>>> 2 eigenvectors failed to convergewarning:eofunc_ts: 2 eigenvectors failed
>>> to convergewarning:avg: Entire input array contains missing values, can't
>>> compute averagewarning:avg: Entire input array contains missing values,
>>> can't compute averagewarning:avg: Entire input array contains missing
>>> values, can't compute averagewarning:avg: Entire input array contains
>>> missing values, can't compute averagewarning:avg: Entire input array
>>> contains missing values, can't compute averagewarning:avg: Entire input
>>> array contains missing values, can't compute averagefatal:Subscript out of
>>> range, error in subscript #0fatal:An error occurred reading
>>> unknownfatal:["Execute.c":8635]:Execute: Error occurred at or near line 207*
>>>
>>> Thank you.
>>>
>>> On Fri, Dec 20, 2019 at 1:15 PM Barry Lynn <barry.h.lynn at gmail.com>
>>> wrote:
>>>
>>>> Hi:
>>>>
>>>> Could you please define what you mean by "missing." Did you
>>>> PrntMinMax(XXXX,False)?
>>>>
>>>>
>>>>
>>>> On Fri, Dec 20, 2019 at 1:20 AM Rahpeni Fajarianti via ncl-talk <
>>>> ncl-talk at ucar.edu> wrote:
>>>>
>>>>> Hi Dennis, Thank you for your time.
>>>>>
>>>>> I am going to get OLR and (850 hPa and 200 hPa) zonal wind anomaly for
>>>>> year 2018 using mjoclivar_2.ncl (
>>>>> https://www.ncl.ucar.edu/Applications/Scripts/mjoclivar_2.ncl) to
>>>>> create netcdf file as input data for mjoclivar_14.
>>>>>
>>>>> The dataset I used were OLR daily mean (
>>>>> https://www.esrl.noaa.gov/psd/data/gridded/data.interp_OLR.html) and
>>>>> u wind daily mean (
>>>>> https://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.pressure.html).
>>>>>
>>>>>
>>>>> When I run the script and I got the result, the wind anomaly data is
>>>>> missing. I used print statement to make sure the values.
>>>>> can anybody help to fix this? Thank you
>>>>>
>>>>> *This is what I got for OLR anomaly:*
>>>>> Variable: f2
>>>>> Type: file
>>>>> filename: olr.day.anomalies.2018
>>>>> path: /home/peni/olr.day.anomalies.2018.nc
>>>>> file global attributes:
>>>>> creation_date : Sel Des 17 08:24:41 WIB 2019
>>>>> source_file : olr.day.mean.nc
>>>>> title : OLR: Daily Anomalies: 2018-2018
>>>>> dimensions:
>>>>> time = 365 // unlimited
>>>>> lat = 73
>>>>> lon = 144
>>>>> variables:
>>>>> double time ( time )
>>>>> units : hours since 1800-01-01 00:00:0.0
>>>>> long_name : Time
>>>>> actual_range : ( 1528872, 1919688 )
>>>>> delta_t : 0000-00-01 00:00:00
>>>>> standard_name : time
>>>>> axis : T
>>>>>
>>>>> float lat ( lat )
>>>>> units : degrees_north
>>>>> actual_range : ( 90, -90 )
>>>>> long_name : Latitude
>>>>> standard_name : latitude
>>>>> axis : Y
>>>>>
>>>>> float lon ( lon )
>>>>> units : degrees_east
>>>>> long_name : Longitude
>>>>> actual_range : ( 0, 360 )
>>>>> standard_name : longitude
>>>>> axis : X
>>>>>
>>>>> float OLR_anom ( time, lat, lon )
>>>>> _FillValue : 32766
>>>>> long_name : Anomalies: Daily OLR
>>>>> units : W/m^2
>>>>>
>>>>> float OLR_anom_sm ( time, lat, lon )
>>>>> _FillValue : 32766
>>>>> long_name : Anomalies from Smooth Daily Climatology
>>>>> units : W/m^2
>>>>>
>>>>>
>>>>> *this is what I got for zonal wind anomaly:*
>>>>> Variable: f3
>>>>> Type: file
>>>>> filename: uwnd.anomalies
>>>>> path: /home/peni/uwnd.anomalies.nc
>>>>> file global attributes:
>>>>> creation_date : Sel Des 17 12:56:50 WIB 2019
>>>>> source_file : uwnd.day.mean.200.nc
>>>>> dimensions:
>>>>> time = 365 // unlimited
>>>>> lat = 73
>>>>> lon = 145
>>>>> variables:
>>>>> double time ( time )
>>>>> long_name : Time
>>>>> units : minutes since 2018-01-01 00:00
>>>>>
>>>>> double lat ( lat )
>>>>> units : degrees_north
>>>>> long_name : Latitude
>>>>>
>>>>> double lon ( lon )
>>>>> units : degrees_east
>>>>> long_name : Longitude
>>>>>
>>>>> float U_anom ( time, lat, lon )
>>>>> long_name : Anomalies from Daily Climatology
>>>>> _FillValue : -9.99e+08
>>>>>
>>>>> float U_anom_sm ( time, lat, lon )
>>>>> long_name : Anomalies from Smooth Daily Climatology
>>>>> _FillValue : -9.99e+08
>>>>>
>>>>>
>>>>> *This is my script:*
>>>>> ;-------------------------------------------------------------
>>>>> ; mjoclivar_2.ncl
>>>>> ;-------------------------------------------------------------
>>>>> ; These files are loaded by default in NCL V6.2.0 and newer
>>>>> ; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
>>>>> ; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
>>>>> ; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
>>>>> ;******************************************************
>>>>>
>>>>> begin
>>>>> ;******************************************************
>>>>> ; User specifications
>>>>> ;******************************************************
>>>>> NC = True ; create netCDF?
>>>>> PLOT = False ; sample plots?
>>>>>
>>>>> ymdStrt = 20180101 ; start yyyymmdd
>>>>> ymdLast = 20181231 ; last
>>>>> yrStrt = ymdStrt/10000
>>>>> yrLast = ymdLast/10000
>>>>>
>>>>> nhar = 4 ; number of fourier comp
>>>>>
>>>>>
>>>>> var = "uwnd2" ; name of file
>>>>>
>>>>> vName = "U" ; name for plots
>>>>>
>>>>> ;diri = "/project/cas/shea/CDC/" ; input directory
>>>>> ;diri = "/Users/shea/Data/AMWG/" ; input directory
>>>>> diri = "/home/peni/" ; new input directory
>>>>> fili = var+".day.mean.850.nc" ; input file
>>>>>
>>>>> if (NC) then
>>>>> diro= "/home/peni/" ; output dir
>>>>> ;filo= var+".day.anomalies.nc" ; output file
>>>>> filo= var+".anomalies.nc" ; output file
>>>>> end if
>>>>>
>>>>> if (PLOT) then
>>>>> wksType = "png" ; send graphics to PNG
>>>>> file
>>>>> wksName = "mjoclivar" ;
>>>>> "mjo."+yrStrt+"_"+yrLast
>>>>> end if
>>>>>
>>>>> ;***********************************************************
>>>>> ; Read user specified time and create required yyyyddd
>>>>>
>>>>> ;***********************************************************
>>>>> f = addfile (diri+fili , "r")
>>>>>
>>>>> time = f->time ; all times on file
>>>>> ymd = cd_calendar(time, -2) ; yyyymmdd
>>>>> iStrt = ind(ymd.eq.ymdStrt) ; index start
>>>>> iLast = ind(ymd.eq.ymdLast) ; index last
>>>>> delete(time)
>>>>> delete(ymd)
>>>>>
>>>>> ;***********************************************************
>>>>> ; Read user specified time and create required yyyyddd
>>>>>
>>>>> ;***********************************************************
>>>>> time = f->time(iStrt:iLast) ; time:units = "hours
>>>>> since"
>>>>> TIME = cd_calendar(time, 0) ; type float
>>>>> year = floattointeger( TIME(:,0) )
>>>>> month = floattointeger( TIME(:,1) )
>>>>> day = floattointeger( TIME(:,2) )
>>>>> ddd = day_of_year(year, month, day)
>>>>> yyyyddd = year*1000 + ddd ; needed for input
>>>>>
>>>>> ;***********************************************************
>>>>> ; Read data: dble2flt
>>>>> ;***********************************************************
>>>>> x = dble2flt( f->$var$(iStrt:iLast,:,:) ) ; convert to
>>>>> float
>>>>> printVarSummary( x )
>>>>>
>>>>> ;***********************************************************
>>>>> ; Compute daily climatology: raw and then 'smoothed'
>>>>> ;***********************************************************
>>>>> xClmDay = clmDayTLL(x, yyyyddd) ; daily climatology at each
>>>>> grid point
>>>>> printVarSummary(xClmDay)
>>>>>
>>>>> ;***********************************************************
>>>>> ; Compute smoothed daily climatology using 'nhar' harmonics
>>>>> ;***********************************************************
>>>>> xClmDay_sm = smthClmDayTLL(xClmDay, nhar)
>>>>> printVarSummary(xClmDay_sm)
>>>>>
>>>>> ;***********************************************************
>>>>> ; Compute daily anomalies using raw and smoothed climatologies
>>>>> ;***********************************************************
>>>>> xAnom = calcDayAnomTLL (x, yyyyddd, xClmDay)
>>>>> printVarSummary(xAnom)
>>>>> printMinMax(xAnom, True)
>>>>>
>>>>> xAnom_sm = calcDayAnomTLL (x, yyyyddd, xClmDay_sm)
>>>>> xAnom_sm at long_name = "Anomalies from Smooth Daily Climatology"
>>>>> printVarSummary(xAnom_sm)
>>>>> printMinMax(xAnom_sm, True)
>>>>>
>>>>> delete( x ) ; no longer needed
>>>>>
>>>>>
>>>>> ;***********************************************************
>>>>> ; Create netCDF: convenience use 'simple' method
>>>>> ;***********************************************************
>>>>>
>>>>> dimx = dimsizes(xAnom)
>>>>> ntim = dimx(0)
>>>>> nlat = dimx(1)
>>>>> mlon = dimx(2)
>>>>>
>>>>> if (NC) then
>>>>>
>>>>> lat = f->lat
>>>>> lon = f->lon
>>>>>
>>>>> system("/bin/rm -f "+diro+filo) ; rm any pre-exist file,
>>>>> if any
>>>>> fnc = addfile (diro+filo, "c")
>>>>>
>>>>> filAtt = 0
>>>>> filAtt at title = vName+": Daily Anomalies: "+yrStrt+"
>>>>> filAtt at source_file = fili
>>>>> filAtt at creation_date = systemfunc("date")
>>>>> fileattdef( fnc, filAtt ) ; copy file attributes
>>>>>
>>>>> setfileoption(fnc,"DefineMode",True)
>>>>>
>>>>> varNC = vName+"_anom"
>>>>> varNC_sm = vName+"_anom_sm"
>>>>>
>>>>> dimNames = (/"time", "lat", "lon"/)
>>>>> dimSizes = (/ -1 , nlat, mlon/)
>>>>> dimUnlim = (/ True , False, False/)
>>>>> filedimdef(fnc,dimNames,dimSizes,dimUnlim)
>>>>>
>>>>> filevardef(fnc, "time" ,typeof(time),getvardims(time))
>>>>> filevardef(fnc, "lat" ,typeof(lat) ,getvardims(lat))
>>>>> filevardef(fnc, "lon" ,typeof(lon) ,getvardims(lon))
>>>>> filevardef(fnc, varNC ,typeof(xAnom) ,getvardims(xAnom))
>>>>>
>>>>> filevardef(fnc, varNC_sm,typeof(xAnom) ,getvardims(xAnom))
>>>>>
>>>>>
>>>>> filevarattdef(fnc,"time" ,time) ; copy
>>>>> time attributes
>>>>> filevarattdef(fnc,"lat" ,lat) ; copy lat
>>>>> attributes
>>>>> filevarattdef(fnc,"lon" ,lon) ; copy lon
>>>>> attributes
>>>>> filevarattdef(fnc,varNC ,xAnom)
>>>>> filevarattdef(fnc,varNC_sm,xAnom_sm)
>>>>>
>>>>> fnc->time = (/time/)
>>>>> fnc->lat = (/lat/)
>>>>> fnc->lon = (/lon/)
>>>>> fnc->$varNC$ = (/xAnom /)
>>>>> fnc->$varNC_sm$ = (/xAnom_sm/)
>>>>> end if
>>>>> end
>>>>>
>>>>>
>>>>>
>>>>> On Thu, Dec 19, 2019 at 9:59 AM Rahpeni Fajarianti <
>>>>> rahpenifajarianti at gmail.com> wrote:
>>>>>
>>>>>>
>>>>>> ---------- Forwarded message ---------
>>>>>> Dari: Rahpeni Fajarianti <rahpenifajarianti at gmail.com>
>>>>>> Date: Kam, 19 Des 2019 9:10 AM
>>>>>> Subject: Fwd: [ncl-talk] fatal:divide: Division by 0, Can't
>>>>>> continuefatal:Div: operator failed, can't continue
>>>>>> fatal:["Execute.c":8635]:Execute: Error occurred at or near line 129
>>>>>> To: reza tisa <rezatisa1 at gmail.com>
>>>>>>
>>>>>>
>>>>>>
>>>>>> ---------- Forwarded message ---------
>>>>>> Dari: Dennis Shea <shea at ucar.edu>
>>>>>> Date: Kam, 19 Des 2019 7:57 AM
>>>>>> Subject: Re: [ncl-talk] fatal:divide: Division by 0, Can't
>>>>>> continuefatal:Div: operator failed, can't continue
>>>>>> fatal:["Execute.c":8635]:Execute: Error occurred at or near line 129
>>>>>> To: Rahpeni Fajarianti <rahpenifajarianti at gmail.com>
>>>>>> Cc: Ncl-talk <ncl-talk at ucar.edu>
>>>>>>
>>>>>>
>>>>>> The GOLDEN RULE of data processing is: ****LOOK AT YOUR DATA****
>>>>>>
>>>>>> This is *user* responsibility. **ALL* the data are 0.0*
>>>>>>
>>>>>> *%>* ncl rdData.ncl
>>>>>> =================
>>>>>> Variable: olr
>>>>>> Dimensions and sizes: [time | 365] x [lat | 13] x [lon | 144]
>>>>>> Coordinates:
>>>>>> time: [1910952..1919688]
>>>>>> lat: [-15..15]
>>>>>> lon: [ 0..357.5]
>>>>>> (0) A*nomalies: Daily OLR (W/m^2) : min=0 max=0*
>>>>>> (0) ------------------------------
>>>>>>
>>>>>> Variable: u200
>>>>>> Dimensions and sizes: [time | 365] x [lat | 13] x [lon | 145]
>>>>>> Coordinates:
>>>>>> time: [ 0..524160]
>>>>>> lat: [ -15.. 15]
>>>>>> lon: [ 0.. 360]
>>>>>> Number Of Attributes: 2
>>>>>> long_name : Anomalies from Daily Climatology
>>>>>> _FillValue : -9.99e+08
>>>>>> (0) *Anomalies from Daily Climatology : min=0 max=0*
>>>>>> (0) ------------------------------
>>>>>>
>>>>>> Variable: u850
>>>>>> Dimensions and sizes: [time | 365] x [lat | 13] x [lon | 145]
>>>>>> Coordinates:
>>>>>> time: [ 0..524160]
>>>>>> lat: [ -15.. 15]
>>>>>> lon: [ 0.. 360]
>>>>>> Number Of Attributes: 2
>>>>>> _FillValue : -9.99e+08
>>>>>>
>>>>>> (0) *Anomalies from Daily Climatology : min=0 max=0*
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> On Wed, Dec 18, 2019 at 4:58 PM Rahpeni Fajarianti <
>>>>>> rahpenifajarianti at gmail.com> wrote:
>>>>>>
>>>>>>> Thank you so much for responding. I attach some files here. The data
>>>>>>> and script that I used.
>>>>>>> The data : olr.day.anomalies.2018.nc , uwnd.anomalies.nc (for
>>>>>>> 200hPa u wind anomalies) , uwnd2.anomalies.nc (for 850hPa u wind
>>>>>>> anomalies)
>>>>>>> my script: mjoclivar_14.ncl
>>>>>>>
>>>>>>> I am very new in NCL and very pleased for help. Thank you.
>>>>>>> mjoclivar_14.ncl
>>>>>>> <https://drive.google.com/file/d/1nstwpHwEVEB9nbuWQPnB9nntN2YVBz79/view?usp=drive_web>
>>>>>>> olr.day.anomalies.2018.nc
>>>>>>> <https://drive.google.com/file/d/1RQIk85UqimEzDzqWZ9V97ZpbqRG8h-cz/view?usp=drive_web>
>>>>>>> uwnd2.anomalies.nc
>>>>>>> <https://drive.google.com/file/d/1ZxF2tZnCr7q1LFQWtGhMOLvaIE3mwK-p/view?usp=drive_web>
>>>>>>> uwnd.anomalies.nc
>>>>>>> <https://drive.google.com/file/d/1rm-f84yy1PkKlA8ZBG87nRGDbnD8auqx/view?usp=drive_web>
>>>>>>>
>>>>>>> On Thu, Dec 19, 2019 at 2:22 AM Dennis Shea <shea at ucar.edu> wrote:
>>>>>>>
>>>>>>>> This is offlineI
>>>>>>>>
>>>>>>>> If you do not know NCL well, please read the excellent tutorial at
>>>>>>>> *http://www.ncl.ucar.edu/Document/Manuals/NCL_User_Guide/*
>>>>>>>> <http://www.ncl.ucar.edu/Document/Manuals/NCL_User_Guide/>
>>>>>>>> ---
>>>>>>>> We [ ncl-talk ] is frustrated by this thread. Please understand,
>>>>>>>> nobody is paid to answer ncl-talk question. We do it on our own time to
>>>>>>>> help people. Debugging is VERY time consuming. In particular, somebody
>>>>>>>> else's code.
>>>>>>>> e
>>>>>>>> The NCL example scripts indicate an approach to use but are not
>>>>>>>> necessarily plug in a file or files and use them. Often, users must make
>>>>>>>> some changes.
>>>>>>>> ---
>>>>>>>> This is the result:
>>>>>>>> Variable: olr
>>>>>>>> Type: float
>>>>>>>> [SNIP]
>>>>>>>> Number Of Attributes: 5
>>>>>>>> _FillValue : 32766
>>>>>>>>
>>>>>>>> (0,0) 32766 ,==== All values are missing [_FillValue ]
>>>>>>>> [SNIP]
>>>>>>>>
>>>>>>>> ===
>>>>>>>> Well .... at what point in your script did things go 'wrong'?
>>>>>>>> Use printVarSummary(), printMinMax(), etc
>>>>>>>> ===
>>>>>>>>
>>>>>>>> I will look at the issue(s) if you make your data available via
>>>>>>>> (say) google drive or ftp.
>>>>>>>>
>>>>>>>> ftp ftp.cgd.ucar.edu
>>>>>>>> anonymous
>>>>>>>> your_email
>>>>>>>> cd incoming
>>>>>>>> put ...
>>>>>>>> put ...
>>>>>>>> put ...
>>>>>>>> quit
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> On Wed, Dec 18, 2019 at 8:14 AM Rahpeni Fajarianti via ncl-talk <
>>>>>>>> ncl-talk at ucar.edu> wrote:
>>>>>>>>
>>>>>>>>> Thanks for your help. I had use print and printVarSummary
>>>>>>>>> statement.
>>>>>>>>>
>>>>>>>>> This is the result:
>>>>>>>>> Variable: olr
>>>>>>>>> Type: float
>>>>>>>>> Total Size: 210240 bytes
>>>>>>>>> 52560 values
>>>>>>>>> Number of Dimensions: 2
>>>>>>>>> Dimensions and sizes: [time | 365] x [lon | 144]
>>>>>>>>> Coordinates:
>>>>>>>>> time: [1910952..1919688]
>>>>>>>>> lon: [ 0..357.5]
>>>>>>>>> Number Of Attributes: 5
>>>>>>>>> _FillValue : 32766
>>>>>>>>> units : W/m^2
>>>>>>>>> long_name : Anomalies: Daily OLR
>>>>>>>>> average_op_ncl : dim_avg_n over dimension(s): lat
>>>>>>>>> wgt_runave_op_ncl : wgt_runave_n
>>>>>>>>> (0,0) 32766
>>>>>>>>> (0,1) 32766
>>>>>>>>> (0,2) 32766
>>>>>>>>> (0,3) 32766
>>>>>>>>> (0,4) 32766
>>>>>>>>> (0,5) 32766
>>>>>>>>> (0,6) 32766
>>>>>>>>> (0,7) 32766
>>>>>>>>> (0,8) 32766
>>>>>>>>> (0,9) 32766
>>>>>>>>> (0,10) 32766
>>>>>>>>> (0,11) 32766
>>>>>>>>> (0,12) 32766
>>>>>>>>> (0,13) 32766
>>>>>>>>> (0,14) 32766
>>>>>>>>> (0,15) 32766
>>>>>>>>> (0,16) 32766
>>>>>>>>> (0,17) 32766
>>>>>>>>> (0,18) 32766
>>>>>>>>> (0,19) 32766
>>>>>>>>> (0,20) 32766
>>>>>>>>> (0,21) 32766
>>>>>>>>> (0,22) 32766
>>>>>>>>> (0,23) 32766
>>>>>>>>> (0,24) 32766
>>>>>>>>> (0,25) 32766
>>>>>>>>> (0,26) 32766
>>>>>>>>> (0,27) 32766
>>>>>>>>> (0,28) 32766
>>>>>>>>> (0,29) 32766
>>>>>>>>> (0,30) 32766
>>>>>>>>> (0,31) 32766
>>>>>>>>> (0,32) 32766
>>>>>>>>> (0,33) 32766
>>>>>>>>> (0,34) 32766
>>>>>>>>> (0,35) 32766
>>>>>>>>> (0,36) 32766
>>>>>>>>> (0,37) 32766
>>>>>>>>> (0,38) 32766
>>>>>>>>> (0,39) 32766
>>>>>>>>> (0,40) 32766
>>>>>>>>> (0,41) 32766
>>>>>>>>> (0,42) 32766
>>>>>>>>> (0,43) 32766
>>>>>>>>> (0,44) 32766
>>>>>>>>> (0,45) 32766
>>>>>>>>> (0,46) 32766
>>>>>>>>> (0,47) 32766
>>>>>>>>> (0,48) 32766
>>>>>>>>> (0,49) 32766
>>>>>>>>> (0,50) 32766
>>>>>>>>> (0,51) 32766
>>>>>>>>> (0,52) 32766
>>>>>>>>> (0,53) 32766
>>>>>>>>> (0,54) 32766
>>>>>>>>> (0,55) 32766
>>>>>>>>> (0,56) 32766
>>>>>>>>> (0,57) 32766
>>>>>>>>> (0,58) 32766
>>>>>>>>> (0,59) 32766
>>>>>>>>> (0,60) 32766
>>>>>>>>> (0,61) 32766
>>>>>>>>> (0,62) 32766
>>>>>>>>> (0,63) 32766
>>>>>>>>> (0,64) 32766
>>>>>>>>> (0,65) 32766
>>>>>>>>> (0,66) 32766
>>>>>>>>> (0,67) 32766
>>>>>>>>> (0,68) 32766
>>>>>>>>> (0,69) 32766
>>>>>>>>> (0,70) 32766
>>>>>>>>> (0,71) 32766
>>>>>>>>> (0,72) 32766
>>>>>>>>> (0,73) 32766
>>>>>>>>> (0,74) 32766
>>>>>>>>> (0,75) 32766
>>>>>>>>> (0,76) 32766
>>>>>>>>> (0,77) 32766
>>>>>>>>> (0,78) 32766
>>>>>>>>> (0,79) 32766
>>>>>>>>> (0,80) 32766
>>>>>>>>> (0,81) 32766
>>>>>>>>> (0,82) 32766
>>>>>>>>> (0,83) 32766
>>>>>>>>> (0,84) 32766
>>>>>>>>> (0,85) 32766
>>>>>>>>> (0,86) 32766
>>>>>>>>> (0,87) 32766
>>>>>>>>> (0,88) 32766
>>>>>>>>> (0,89) 32766
>>>>>>>>> (0,90) 32766
>>>>>>>>> (0,91) 32766
>>>>>>>>> (0,92) 32766
>>>>>>>>> (0,93) 32766
>>>>>>>>> (0,94) 32766
>>>>>>>>> (0,95) 32766
>>>>>>>>> (0,96) 32766
>>>>>>>>> (0,97) 32766
>>>>>>>>> (0,98) 32766
>>>>>>>>> (0,99) 32766
>>>>>>>>> (0,100) 32766
>>>>>>>>> (0,101) 32766
>>>>>>>>> (0,102) 32766
>>>>>>>>> (0,103) 32766
>>>>>>>>> (0,104) 32766
>>>>>>>>> (0,105) 32766
>>>>>>>>> (0,106) 32766
>>>>>>>>> (0,107) 32766
>>>>>>>>> (0,108) 32766
>>>>>>>>> (0,109) 32766
>>>>>>>>> (0,110) 32766
>>>>>>>>> (0,111) 32766
>>>>>>>>> (0,112) 32766
>>>>>>>>> (0,113) 32766
>>>>>>>>> (0,114) 32766
>>>>>>>>> (0,115) 32766
>>>>>>>>> (0,116) 32766
>>>>>>>>> (0,117) 32766
>>>>>>>>> (0,118) 32766
>>>>>>>>> (0,119) 32766
>>>>>>>>> (0,120) 32766
>>>>>>>>> (0,121) 32766
>>>>>>>>> (0,122) 32766
>>>>>>>>> (0,123) 32766
>>>>>>>>> (0,124) 32766
>>>>>>>>> (0,125) 32766
>>>>>>>>> (0,126) 32766
>>>>>>>>> (0,127) 32766
>>>>>>>>> (0,128) 32766
>>>>>>>>> (0,129) 32766
>>>>>>>>> (0,130) 32766
>>>>>>>>> (0,131) 32766
>>>>>>>>> (0,132) 32766
>>>>>>>>> (0,133) 32766
>>>>>>>>> (0,134) 32766
>>>>>>>>> (0,135) 32766
>>>>>>>>> (0,136) 32766
>>>>>>>>> (0,137) 32766
>>>>>>>>> (0,138) 32766 ect..
>>>>>>>>>
>>>>>>>>> What should I fix to my script?
>>>>>>>>> Thanks.
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> On Wed, Dec 18, 2019 at 8:50 AM Adam Phillips <asphilli at ucar.edu>
>>>>>>>>> wrote:
>>>>>>>>>
>>>>>>>>>> Hi Rahpeni,
>>>>>>>>>> Please always include ncl-talk on your replies; that way others
>>>>>>>>>> can assist and see the correspondence.
>>>>>>>>>>
>>>>>>>>>> You stated:
>>>>>>>>>> Thank you so much, it helps !
>>>>>>>>>>
>>>>>>>>>> I am still getting some errors.
>>>>>>>>>> I got :
>>>>>>>>>> fatal:Subscript out of range, error in subscript #1
>>>>>>>>>> fatal:An error occurred reading olr
>>>>>>>>>> fatal:["Execute.c":8635]:Execute: Error occurred at or near line
>>>>>>>>>> 141
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> Note that it is for you to attempt to diagnose the issue before
>>>>>>>>>> responding to ncl-talk. The error message is telling you what the problem
>>>>>>>>>> is: You are attempting to access an index of a dimension that is bigger
>>>>>>>>>> than the dimension itself. The lines around 141 are (from the email you
>>>>>>>>>> sent me):
>>>>>>>>>>
>>>>>>>>>> ncl 140> do ml=0,mlon-1
>>>>>>>>>> ncl 141> cdata(ml ,:) = (/ olr(:,ml) /)
>>>>>>>>>> ncl 142> cdata(ml+ mlon,:) = (/ u850(:,ml) /)
>>>>>>>>>> ncl 143> cdata(ml+2*mlon,:) = (/ u200(:,ml) /)
>>>>>>>>>> ncl 144> end do
>>>>>>>>>>
>>>>>>>>>> NCL is zero based, so error in subscript #1 is referring to the
>>>>>>>>>> 2nd index, line 141 shows that you are accessing the 2nd index of olr like
>>>>>>>>>> this: (/ olr(:,ml) /)
>>>>>>>>>>
>>>>>>>>>> ml at some point is bigger than the 2nd dimension size of olr.
>>>>>>>>>> You will have to figure out why. Use print and printVarSummary statements
>>>>>>>>>> to help to diagnose the issue.
>>>>>>>>>>
>>>>>>>>>> Hope that helps!
>>>>>>>>>> Adam
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> On Tue, Dec 17, 2019 at 1:40 AM Rahpeni Fajarianti via ncl-talk <
>>>>>>>>>> ncl-talk at ucar.edu> wrote:
>>>>>>>>>>
>>>>>>>>>>> Hi all !
>>>>>>>>>>> I want to make MJO clivar 14, but I have some problem.
>>>>>>>>>>> This is :
>>>>>>>>>>> fatal:divide: Division by 0, Can't continue
>>>>>>>>>>> fatal:Div: operator failed, can't continue
>>>>>>>>>>> fatal:["Execute.c":8635]:Execute: Error occurred at or near line
>>>>>>>>>>> 129
>>>>>>>>>>>
>>>>>>>>>>> This is my script:
>>>>>>>>>>> ncl 0> load
>>>>>>>>>>> "$NCARG_ROOT/lib/ncarg/nclscripts/csm/diagnostics_cam.ncl"
>>>>>>>>>>> ncl 1> ;******************************************************
>>>>>>>>>>> u850 = dim_rmvmean_n(u850, 0)
>>>>>>>>>>> u200 = dim_rmvmean_nncl 2> ;
>>>>>>>>>>> (uncl 3> ; mjoclivar_14.ncl
>>>>>>>>>>> ncl 4> ;
>>>>>>>>>>> ncl 5>
>>>>>>>>>>> ;***********************************************************
>>>>>>>>>>> ncl 6> ; Combined EOFs
>>>>>>>>>>> hncl 7> ; Latest Update: July, 2016: Eun-Pa Lim; Bureau of
>>>>>>>>>>> Meteorology, Australia
>>>>>>>>>>> ncl 8>
>>>>>>>>>>> ;***********************************************************
>>>>>>>>>>> iancncl 9> ;;
>>>>>>>>>>> nncl 10> ;; The following are automatically loaded from
>>>>>>>>>>> 6.2.0 onward
>>>>>>>>>>> ncl 11> ;;load
>>>>>>>>>>> "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
>>>>>>>>>>> ncl 12> ;;load
>>>>>>>>>>> "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
>>>>>>>>>>> ncl 13> ;;load
>>>>>>>>>>> "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
>>>>>>>>>>> ncl 14>
>>>>>>>>>>> ncl 15> undef("read_rename")
>>>>>>>>>>> ncl 16> function read_rename(f[1]:file, varName[1]:string \
>>>>>>>>>>> ncl 16> ,iStrt[1]:integer, iLast[1]:integer \
>>>>>>>>>>> ncl 16> ,latS[1]:numeric , latN[1]:numeric )
>>>>>>>>>>> ncl 17> ; Utility to force specific named dimensions
>>>>>>>>>>> ncl 18> ; This is done for historical reasons (convenience)
>>>>>>>>>>> ncl 19> begin
>>>>>>>>>>> ncl 20> work = f->$varName$(iStrt:iLast,{latS:latN},:) ;
>>>>>>>>>>> (time,lat,lon)
>>>>>>>>>>> ncl 21> work!0 = "time" ;
>>>>>>>>>>> CAM model names
>>>>>>>>>>> ncl 22> work!1 = "lat"
>>>>>>>>>>> ncl 23> work!2 = "lon"
>>>>>>>>>>> ncl 24> return(work)
>>>>>>>>>>> ncl 25> end
>>>>>>>>>>> ncl 26> ; =========================> MAIN
>>>>>>>>>>> <==============================
>>>>>>>>>>> ncl 27> begin
>>>>>>>>>>> ncl 28> neof = 2
>>>>>>>>>>> ncl 29>
>>>>>>>>>>> ncl 30> latS = -15
>>>>>>>>>>> ncl 31> latN = 15
>>>>>>>>>>> ncl 32>
>>>>>>>>>>> ncl 33> ymdStrt = 20180101 ; start
>>>>>>>>>>> yyyymmdd
>>>>>>>>>>> ncl 34> ymdLast = 20181231 ; last
>>>>>>>>>>> ncl 35>
>>>>>>>>>>> ncl 36> yrStrt = ymdStrt/10000
>>>>>>>>>>> ncl 37> yrLast = ymdLast/10000
>>>>>>>>>>> ncl 38>
>>>>>>>>>>> ncl 39> pltDir = "/home/peni/" ;
>>>>>>>>>>> plot directory
>>>>>>>>>>> *ncl 40> pltType = "png"
>>>>>>>>>>> ncl 41> pltName = "mjoclivar"
>>>>>>>>>>> ncl 42>
>>>>>>>>>>> ,ncl 43> diri = "/home/peni/"
>>>>>>>>>>> ; input directory
>>>>>>>>>>> ncl 44>
>>>>>>>>>>> ncl 45> filolr = "olr.day.anomalies.2018.nc"
>>>>>>>>>>> ncl 46> filu200 = "uwnd.anomalies.nc"
>>>>>>>>>>> ncl 47> filu850 = "uwnd2.anomalies.nc.nc"
>>>>>>>>>>> ncl 48>
>>>>>>>>>>> ncl 49> ;************************************************
>>>>>>>>>>> ncl 50> ; create BandPass Filter
>>>>>>>>>>> ncl 51> ;************************************************
>>>>>>>>>>> ncl 52> ihp = 2 ; bpf=>band
>>>>>>>>>>> pass filter
>>>>>>>>>>> ncl 53> nWgt = 201
>>>>>>>>>>> ncl 54> sigma = 1.0 ; Lanczos
>>>>>>>>>>> sigma
>>>>>>>>>>> ncl 55> fca = 1./100.
>>>>>>>>>>> ncl 56> fcb = 1./20.
>>>>>>>>>>> ncl 57> wgt = filwgts_lanczos (nWgt, ihp, fca, fcb, sigma
>>>>>>>>>>> )
>>>>>>>>>>> ncl 58>
>>>>>>>>>>> )ncl 59>
>>>>>>>>>>> ;***********************************************************
>>>>>>>>>>> ncl 60> ; Find the indices corresponding to the start/end times
>>>>>>>>>>> ncl 61>
>>>>>>>>>>> ;***********************************************************
>>>>>>>>>>> ncl 62> f = addfile (diri+filolr , "r")
>>>>>>>>>>>
>>>>>>>>>>> ncl 63> TIME = f->time ; days
>>>>>>>>>>> since ...
>>>>>>>>>>> ncl 64>
>>>>>>>>>>> ncl 65> YMD = cd_calendar(TIME, -2) ; entire
>>>>>>>>>>> (time,6)
>>>>>>>>>>> ncl 66>
>>>>>>>>>>> ncl 67> iStrt = ind(YMD.eq.ymdStrt) ; index
>>>>>>>>>>> start
>>>>>>>>>>> ncl 68> iLast = ind(YMD.eq.ymdLast) ; index
>>>>>>>>>>> last
>>>>>>>>>>> ncl 69> delete([/ TIME, YMD /])
>>>>>>>>>>> ncl 70>
>>>>>>>>>>> ncl 71>
>>>>>>>>>>> ;***********************************************************
>>>>>>>>>>> ncl 72> ; Read anomalies
>>>>>>>>>>> ncl 73>
>>>>>>>>>>> ;***********************************************************
>>>>>>>>>>> ncl 74>
>>>>>>>>>>> ncl 75> work =
>>>>>>>>>>> read_rename(f,"OLR_anom",iStrt,iLast,latS,latN) ; (time,lat,lon)
>>>>>>>>>>> ncl 76> OLR = dim_avg_n_Wrap(work, 1)
>>>>>>>>>>> ; (time,lon)
>>>>>>>>>>> ncl 77> delete(work)
>>>>>>>>>>> ncl 78>
>>>>>>>>>>> ncl 79> f = addfile (diri+filu850 , "r")
>>>>>>>>>>>
>>>>>>>>>>> ncl 80> work =
>>>>>>>>>>> read_rename(f,"U_anom",iStrt,iLast,latS,latN) ; (time,lat,lon)
>>>>>>>>>>> ncl 81> U850 = dim_avg_n_Wrap(work, 1) ;
>>>>>>>>>>> (time,lon)
>>>>>>>>>>> ncl 82> delete(work)
>>>>>>>>>>> ncl 83>
>>>>>>>>>>> ncl 84> f = addfile (diri+filu200 , "r")
>>>>>>>>>>>
>>>>>>>>>>> ncl 85> work =
>>>>>>>>>>> read_rename(f,"U_anom",iStrt,iLast,latS,latN) ; (time,lat,lon)
>>>>>>>>>>> ncl 86> U200 = dim_avg_n_Wrap(work, 1) ;
>>>>>>>>>>> (time,lon)
>>>>>>>>>>> ncl 87>
>>>>>>>>>>> ncl 88> dimw = dimsizes( work )
>>>>>>>>>>> ncl 89> ntim = dimw(0)
>>>>>>>>>>> ncl 90> nlat = dimw(1)
>>>>>>>>>>> ncl 91> mlon = dimw(2)
>>>>>>>>>>> ncl 92> delete(work)
>>>>>>>>>>> ncl 93>
>>>>>>>>>>> ncl 94> lon = OLR&lon
>>>>>>>>>>>
>>>>>>>>>>> ncl 95> time = OLR&time
>>>>>>>>>>> ncl 96> date = cd_calendar(time, -2) ; yyyymmdd
>>>>>>>>>>> ncl 97>
>>>>>>>>>>> ncl 98> ;************************************************
>>>>>>>>>>> ncl 99> ; Apply the band pass filter to the original anomalies
>>>>>>>>>>> ncl 100> ;************************************************
>>>>>>>>>>> ncl 101> olr = wgt_runave_n_Wrap ( OLR, wgt, 0, 0) ;
>>>>>>>>>>> (time,lon)
>>>>>>>>>>> )
>>>>>>>>>>> ;************************************************
>>>>>>>>>>>
>>>>>>>>>>> imax_olr_eof1 = maxind(ceof(0,0,:))
>>>>>>>>>>> imax_olr_eof2 ncl 102> u850 = wgt_runave_n_Wrap (U850,
>>>>>>>>>>> wgt, 0, 0)
>>>>>>>>>>> ncl 103> u200 = wgt_runave_n_Wrap (U200, wgt, 0, 0)
>>>>>>>>>>> ncl 104>
>>>>>>>>>>> ncl 105> ;************************************************
>>>>>>>>>>> ncl 106> ; remove temporal means of band pass series: *not*
>>>>>>>>>>> necessary
>>>>>>>>>>> ncl 107> ;************************************************
>>>>>>>>>>> ncl 108> olr = dim_rmvmean_n( olr, 0) ;
>>>>>>>>>>> (time,lon)
>>>>>>>>>>> ncl 109> u850 = dim_rmvmean_n(u850, 0)
>>>>>>>>>>> ncl 110> u200 = dim_rmvmean_n(u200, 0)
>>>>>>>>>>> ncl 111>
>>>>>>>>>>> ncl 112> ;************************************************
>>>>>>>>>>> ncl 113> ; Compute the temporal variance at each lon
>>>>>>>>>>> ncl 114> ;************************************************
>>>>>>>>>>> ncl 115> var_olr = dim_variance_n_Wrap( olr, 0) ; (lon)
>>>>>>>>>>> ncl 116> var_u850 = dim_variance_n_Wrap(u850, 0)
>>>>>>>>>>> ncl 117> var_u200 = dim_variance_n_Wrap(u200, 0)
>>>>>>>>>>> ncl 118>
>>>>>>>>>>> ncl 119> ;************************************************
>>>>>>>>>>> ncl 120> ; Compute the zonal mean of the temporal variance
>>>>>>>>>>> Oncl 121> ;************************************************
>>>>>>>>>>> ncl 122> zavg_var_olr = dim_avg_n_Wrap( var_olr , 0)
>>>>>>>>>>> ncl 123> zavg_var_u850 = dim_avg_n_Wrap( var_u850, 0)
>>>>>>>>>>> ncl 124> zavg_var_u200 = dim_avg_n_Wrap( var_u200, 0)
>>>>>>>>>>> ncl 125>
>>>>>>>>>>> ncl 126> ;************************************************
>>>>>>>>>>> ncl 127> ; Normalize by sqrt(avg_var*)
>>>>>>>>>>> ncl 128> ;************************************************
>>>>>>>>>>> ncl 129> olr = olr/sqrt(zavg_var_olr ) ; (time,lon)
>>>>>>>>>>> ncl 130> u850 = u850/sqrt(zavg_var_u850)
>>>>>>>>>>> ncl 131> u200 = u200/sqrt(zavg_var_u200)
>>>>>>>>>>> ncl 132>
>>>>>>>>>>> cncl 133> ;************************************************
>>>>>>>>>>> ncl 134> ; Combine the normalized data into one variable
>>>>>>>>>>> ncl 135> ;************************************************
>>>>>>>>>>> ncl 136> cdata = new ( (/3*mlon,ntim/), typeof(olr),
>>>>>>>>>>> getFillValue(olr))
>>>>>>>>>>> ncl 137> do ml=0,mlon-1
>>>>>>>>>>> ncl 138> cdata(ml ,:) = (/ olr(:,ml) /)
>>>>>>>>>>> ncl 139> cdata(ml+ mlon,:) = (/ u850(:,ml) /)
>>>>>>>>>>> ncl 140> cdata(ml+2*mlon,:) = (/ u200(:,ml) /)
>>>>>>>>>>> ncl 141> end do
>>>>>>>>>>> ncl 142>
>>>>>>>>>>> =ncl 143> ;************************************************
>>>>>>>>>>> ncl 144> ; Compute **combined** EOF; Sign of EOF is arbitrary
>>>>>>>>>>> ncl 145> ;************************************************
>>>>>>>>>>> ncl 146> eof_cdata = eofunc(cdata , neof, False) ;
>>>>>>>>>>> (neof,3*mlon)
>>>>>>>>>>> ncl 147> print("==============")
>>>>>>>>>>> ncl 148> printVarSummary(eof_cdata)
>>>>>>>>>>> ncl 149> printMinMax(eof_cdata, True)
>>>>>>>>>>> ncl 150>
>>>>>>>>>>> ncl 151> eof_ts_cdata = eofunc_ts(cdata,eof_cdata,False) ;
>>>>>>>>>>> (neof,3*ntim)
>>>>>>>>>>> ncl 152> print("==============")
>>>>>>>>>>>
>>>>>>>>>>> ncl 153> printVarSummary(eof_ts_cdata)
>>>>>>>>>>> ncl 154> printMinMax(eof_ts_cdata, True)
>>>>>>>>>>> ncl 155>
>>>>>>>>>>> incl 156> ;************************************************
>>>>>>>>>>> ncl 157> ; For clarity, explicitly extract each variable. Create
>>>>>>>>>>> time series
>>>>>>>>>>> ncl 158> ;************************************************
>>>>>>>>>>> ncl 159>
>>>>>>>>>>> tncl 160> nvar = 3 ; "olr", "u850", "u200"
>>>>>>>>>>> ncl 161> ceof = new( (/nvar,neof,mlon/), typeof(cdata),
>>>>>>>>>>> getFillValue(cdata))
>>>>>>>>>>> ncl 162>
>>>>>>>>>>> dncl 163> do n=0,neof-1
>>>>>>>>>>> ncl 164> ceof(0,n,:) = eof_cdata(n,0:mlon-1) ; olr
>>>>>>>>>>> ncl 165> ceof(1,n,:) = eof_cdata(n,mlon:2*mlon-1) ; u850
>>>>>>>>>>> ncl 166> ceof(2,n,:) = eof_cdata(n,2*mlon:) ; u200
>>>>>>>>>>> ncl 167> end do
>>>>>>>>>>> ncl 168>
>>>>>>>>>>> (ncl 169> ceof!0 = "var"
>>>>>>>>>>> ncl 170> ceof!1 = "eof"
>>>>>>>>>>> ncl 171> ceof!2 = "lon"
>>>>>>>>>>> ncl 172> ceof&lon = olr&lon
>>>>>>>>>>> ncl 173>
>>>>>>>>>>> incl 174> ceof_ts = new( (/nvar,neof,ntim/),
>>>>>>>>>>> typeof(cdata), getFillValue(cdata))
>>>>>>>>>>> ncl 175> ceof_ts(0,:,:) = eofunc_ts_Wrap(
>>>>>>>>>>> olr(lon|:,time|:),ceof(0,:,:),False) ; (0,neof,ntim)
>>>>>>>>>>> ncl 176> ceof_ts(1,:,:) =
>>>>>>>>>>> eofunc_ts_Wrap(u850(lon|:,time|:),ceof(1,:,:),False) ; (1,neof,ntim)
>>>>>>>>>>> ncl 177> ceof_ts(2,:,:) =
>>>>>>>>>>> eofunc_ts_Wrap(u200(lon|:,time|:),ceof(2,:,:),False) ; (2,neof,ntim)
>>>>>>>>>>> ncl 178>
>>>>>>>>>>> ncl 179> ;**********************************************t*
>>>>>>>>>>> ncl 180> ; Add code contributed by Marcus N. Morgan, Florida
>>>>>>>>>>> Institute of Technology; Feb 2015
>>>>>>>>>>> ncl 181> ; Calculate % variance (pcv_ )accounted for by OLR,
>>>>>>>>>>> U850 and U200
>>>>>>>>>>> ncl 182> ;************************************************
>>>>>>>>>>> ncl 183>
>>>>>>>>>>> ncl 184> pcv_eof_olr = new(neof,typeof(ceof))
>>>>>>>>>>> ncl 185> pcv_eof_u850 = new(neof,typeof(ceof))
>>>>>>>>>>> ncl 186> pcv_eof_u200 = new(neof,typeof(ceof))
>>>>>>>>>>> ncl 187>
>>>>>>>>>>> ncl 188> do n=0,neof-1
>>>>>>>>>>> ncl 189> pcv_eof_olr(n) = avg((ceof(0,n,:)*sqrt(ceof at eval
>>>>>>>>>>> (n)))^2)*100
>>>>>>>>>>> ncl 190> pcv_eof_u850(n) = avg((ceof(1,n,:)*sqrt(ceof at eval
>>>>>>>>>>> (n)))^2)*100
>>>>>>>>>>> ncl 191> pcv_eof_u200(n) = avg((ceof(2,n,:)*sqrt(ceof at eval
>>>>>>>>>>> (n)))^2)*100
>>>>>>>>>>> ncl 192> ;;print("pcv: neof="+(n+1)+": "+pcv_eof_olr(n)+"
>>>>>>>>>>> "+pcv_eof_u850(n)+" "+pcv_eof_u200(n))
>>>>>>>>>>> ncl 193> end do
>>>>>>>>>>> ncl 194>
>>>>>>>>>>> ncl 195> ;************************************************
>>>>>>>>>>> ncl 196> ; Change sign of EOFs for spatial structures of PC1 and
>>>>>>>>>>> PC2
>>>>>>>>>>> ncl 197> ; to represent convection over the tropical Indian
>>>>>>>>>>> Ocean and the tropical western Pacific Ocean, respectively
>>>>>>>>>>> ncl 198> ; (Ad hoc approach)
>>>>>>>>>>> ncl 199> ;************************************************
>>>>>>>>>>> ncl 200>
>>>>>>>>>>> ncl 201> imax_olr_eof1 = maxind(ceof(0,0,:))
>>>>>>>>>>> ncl 202> imax_olr_eof2 = maxind(ceof(0,1,:))
>>>>>>>>>>> ncl 203>
>>>>>>>>>>> )ncl 204> lonmax_eof1 = ceof&lon(imax_olr_eof1) ;
>>>>>>>>>>> longitude of max value (i.e. suppressed convection)
>>>>>>>>>>> ncl 205> lonmax_eof2 = ceof&lon(imax_olr_eof2)
>>>>>>>>>>> ncl 206>
>>>>>>>>>>> sncl 207> if (lonmax_eof1.ge.60 .and. lonmax_eof1.lt.180) then
>>>>>>>>>>> ; Change the sign of EOF1
>>>>>>>>>>> ncl 208> ceof(:,0,:) = -ceof(:,0,:)
>>>>>>>>>>> ; if OLR is positive
>>>>>>>>>>> ncl 209> ceof_ts(:,0,:) = -ceof_ts(:,0,:)
>>>>>>>>>>> ; over the tropical Indian Ocean
>>>>>>>>>>> ncl 210> eof_cdata(0,:) = -eof_cdata(0,:)
>>>>>>>>>>> ncl 211> eof_ts_cdata(0,:) = -eof_ts_cdata(0,:)
>>>>>>>>>>> ncl 212> end if
>>>>>>>>>>> ncl 213>
>>>>>>>>>>> encl 214> if (lonmax_eof2.ge.120 .and. lonmax_eof2.lt.180)
>>>>>>>>>>> then ; Change the sign of EOF2
>>>>>>>>>>> ncl 215> ceof(:,1,:) = -ceof(:,1,:)
>>>>>>>>>>> ; if OLR is positive
>>>>>>>>>>> ncl 216> ceof_ts(:,1,:) = -ceof_ts(:,1,:)
>>>>>>>>>>> ; over the tropical western Pacific Ocean
>>>>>>>>>>> ncl 217> eof_cdata(1,:) = -eof_cdata(1,:)
>>>>>>>>>>> ncl 218> eof_ts_cdata(1,:) = -eof_ts_cdata(1,:)
>>>>>>>>>>> ncl 219> end if
>>>>>>>>>>> ncl 220>
>>>>>>>>>>> ncl 221> print("==============")
>>>>>>>>>>> rncl 222> printVarSummary(eof_cdata)
>>>>>>>>>>> ncl 223> printMinMax(eof_cdata, True)
>>>>>>>>>>> ncl 224>
>>>>>>>>>>> ncl 225> ;************************************************
>>>>>>>>>>> ncl 226> ; Compute cross correlation of each variable's EOF time
>>>>>>>>>>> series at zero-lag
>>>>>>>>>>> ncl 227> ;************************************************
>>>>>>>>>>> ncl 228> r_olr_u850 = escorc(ceof_ts(0,:,:) , ceof_ts(1,:,:)
>>>>>>>>>>> ) ; (neof)
>>>>>>>>>>> ncl 229> r_olr_u200 = escorc(ceof_ts(0,:,:) , ceof_ts(2,:,:) )
>>>>>>>>>>> ncl 230> r_u850_u200 = escorc(ceof_ts(1,:,:) , ceof_ts(2,:,:) )
>>>>>>>>>>> ncl 231>
>>>>>>>>>>> sncl 232> print("==============")
>>>>>>>>>>> ncl 233> do n=0,neof-1
>>>>>>>>>>> ncl 234> print("neof="+n \
>>>>>>>>>>> ncl 234> +" r_olr_u850="
>>>>>>>>>>> +sprintf("%4.3f",r_olr_u850(n)) \
>>>>>>>>>>> ncl 234> +" r_olr_u200="
>>>>>>>>>>> +sprintf("%4.3f",r_olr_u200(n)) \
>>>>>>>>>>> ncl 234> +"
>>>>>>>>>>> r_u850_u200="+sprintf("%4.3f",r_u850_u200(n)) )
>>>>>>>>>>> ncl 235> end do
>>>>>>>>>>> ncl 236> print("==============")
>>>>>>>>>>> ncl 237>
>>>>>>>>>>> ncl 238> ;************************************************
>>>>>>>>>>> 5N: "+yrStrt+"-"+yrLast
>>>>>>>>>>>
>>>>>>>>>>> do ncl 239> ; Compute cross correlation of the multivariate
>>>>>>>>>>> EOF; EOF 1 vs EOF 2
>>>>>>>>>>> ncl 240> ;************************************************
>>>>>>>>>>> ncl 241>
>>>>>>>>>>> ncl 242> mxlag = 25
>>>>>>>>>>> ncl 243> rlag_01 =
>>>>>>>>>>> esccr(eof_ts_cdata(0,:),eof_ts_cdata(1,:), mxlag) ; (N,mxlag+1)
>>>>>>>>>>> ncl 244> rlag_10 =
>>>>>>>>>>> esccr(eof_ts_cdata(1,:),eof_ts_cdata(0,:), mxlag) ; (N,mxlag+1)
>>>>>>>>>>> ncl 245> ccr_12 = new ( (/2*mxlag+1/), float)
>>>>>>>>>>> ncl 246>
>>>>>>>>>>> cncl 247> ccr_12(mxlag:) = rlag_10(0:mxlag)
>>>>>>>>>>> ncl 248> ccr_12(0:mxlag) = rlag_01(::-1) ; reverse
>>>>>>>>>>> order
>>>>>>>>>>> ncl 249> ;;print(ccr_12)
>>>>>>>>>>> ncl 250>
>>>>>>>>>>> sncl 251>
>>>>>>>>>>> Pncl 252> ;************************************************
>>>>>>>>>>> ncl 253> ; Normalize the multivariate EOF 1&2 component time
>>>>>>>>>>> series
>>>>>>>>>>> ncl 254> ; Compute (PC1^2+PC2^2): values > 1 indicate "strong"
>>>>>>>>>>> periods
>>>>>>>>>>> ncl 255> ;************************************************
>>>>>>>>>>> ncl 256> eof_ts_cdata(0,:) =
>>>>>>>>>>> eof_ts_cdata(0,:)/stddev(eof_ts_cdata(0,:))
>>>>>>>>>>> ncl 257> eof_ts_cdata(1,:) =
>>>>>>>>>>> eof_ts_cdata(1,:)/stddev(eof_ts_cdata(1,:))
>>>>>>>>>>> ncl 258>
>>>>>>>>>>> ncl 259> mjo_ts_index = eof_ts_cdata(0,:)^2 +
>>>>>>>>>>> eof_ts_cdata(1,:)^2
>>>>>>>>>>> ncl 260> mjo_ts_index_smt = runave(mjo_ts_index, 91, 0) ;
>>>>>>>>>>> 91-day running mean
>>>>>>>>>>> ncl 261>
>>>>>>>>>>> ncl 262> nGood = num(.not.ismissing(mjo_ts_index)) ; #
>>>>>>>>>>> non-missing
>>>>>>>>>>> ncl 263> nStrong = num(mjo_ts_index .ge. 1.0)
>>>>>>>>>>> ncl 264> print("nGood="+nGood+" nStrong="+nStrong+"
>>>>>>>>>>> nOther="+(nGood-nStrong))
>>>>>>>>>>> ncl 265>
>>>>>>>>>>> rncl 266> ;************************************************
>>>>>>>>>>> ncl 267> ; Write PC results to netCDF for use in another example.
>>>>>>>>>>> ncl 268> ;************************************************
>>>>>>>>>>> ncl 269> mjo_ts_index!0 = "time"
>>>>>>>>>>> ncl 270> mjo_ts_index&time = time
>>>>>>>>>>> ncl 271> mjo_ts_index at long_name = "MJO PC INDEX"
>>>>>>>>>>> ncl 272> mjo_ts_index at info = "(PC1^2 + PC2^2)"
>>>>>>>>>>> ncl 273>
>>>>>>>>>>> ncl 274> PC1 = eof_ts_cdata(0,:)
>>>>>>>>>>> ncl 275> PC1!0 = "time"
>>>>>>>>>>> ncl 276> PC1&time = time
>>>>>>>>>>> ncl 277> PC1 at long_name = "PC1"
>>>>>>>>>>> ncl 278> PC1 at info = "PC1/stddev(PC1)"
>>>>>>>>>>> ncl 279>
>>>>>>>>>>> oncl 280> PC2 = eof_ts_cdata(1,:)
>>>>>>>>>>> ncl 281> PC2!0 = "time"
>>>>>>>>>>> ncl 282> PC2&time = time
>>>>>>>>>>> ncl 283> PC2 at long_name = "PC2"
>>>>>>>>>>> ncl 284> PC2 at info = "PC2/stddev(PC2)"
>>>>>>>>>>> ncl 285>
>>>>>>>>>>> ncl 286> diro = "./"
>>>>>>>>>>> ncl 287> filo = "MJO_PC_INDEX.nc"
>>>>>>>>>>> ncl 288> system("/bin/rm -f "+diro+filo) ; remove any
>>>>>>>>>>> pre-existing file
>>>>>>>>>>> ncl 289> ncdf = addfile(diro+filo,"c") ; open output
>>>>>>>>>>> netCDF file
>>>>>>>>>>> of legend.
>>>>>>>>>>> ncl 290> ; make time an
>>>>>>>>>>> UNLIMITED dimension
>>>>>>>>>>> ncl 291> filedimdef(ncdf,"time",-1,True) ; recommended for
>>>>>>>>>>> most applications
>>>>>>>>>>> ncl 292> ; output variables
>>>>>>>>>>> directly
>>>>>>>>>>> ncl 293> ncdf->MJO_INDEX = mjo_ts_index
>>>>>>>>>>> ncl 294> ncdf->PC1 = PC1
>>>>>>>>>>> ncl 295> ncdf->PC2 = PC2
>>>>>>>>>>> ncl 296>
>>>>>>>>>>> ncl 297>
>>>>>>>>>>> ;------------------------------------------------------------
>>>>>>>>>>> ncl 298> ; PLOTS
>>>>>>>>>>> ncl 299>
>>>>>>>>>>> ;------------------------------------------------------------
>>>>>>>>>>> ncl 300>
>>>>>>>>>>> ncl 301> yyyymmdd = cd_calendar(time, -2)
>>>>>>>>>>> tncl 302> yrfrac = yyyymmdd_to_yyyyfrac(yyyymmdd, 0.0)
>>>>>>>>>>> ncl 303> delete([/ yrfrac at long_name, lon at long_name /])
>>>>>>>>>>> ncl 304>
>>>>>>>>>>> ncl 305> day = ispan(-mxlag, mxlag, 1)
>>>>>>>>>>> ncl 306> ;day at long_name = "lag (day)"
>>>>>>>>>>> ncl 307>
>>>>>>>>>>> ncl 308> pltPath = pltDir+pltName
>>>>>>>>>>> ncl 309>
>>>>>>>>>>> ncl 310> wks = gsn_open_wks(pltType,pltPath)
>>>>>>>>>>> ncl 311> gsn_define_colormap(wks,"default")
>>>>>>>>>>> ncl 312> plot = new(3,graphic)
>>>>>>>>>>> ncl 313>
>>>>>>>>>>> ncl 314> ;************************************************
>>>>>>>>>>> ncl 315> ; Multivariate EOF plots
>>>>>>>>>>> ncl 316> ;************************************************
>>>>>>>>>>> ncl 317> rts = True
>>>>>>>>>>> ncl 318> rts at gsnDraw = False ; don't draw yet
>>>>>>>>>>> ncl 319> rts at gsnFrame = False ; don't advance frame yet
>>>>>>>>>>> ncl 320> rts at gsnScale = True ; force text scaling
>>>>>>>>>>>
>>>>>>>>>>> ncl 321>
>>>>>>>>>>> ncl 322> rts at vpHeightF = 0.40 ; Changes the aspect
>>>>>>>>>>> ratio
>>>>>>>>>>> ncl 323> rts at vpWidthF = 0.85
>>>>>>>>>>> ncl 324> rts at vpXF = 0.10 ; change start locations
>>>>>>>>>>> ncl 325> rts at vpYF = 0.75 ; the plot
>>>>>>>>>>> ncl 326> rts at xyLineThicknesses = (/2, 2, 2/)
>>>>>>>>>>> ncl 327> rts at xyLineColors = (/"black","red","blue"/)
>>>>>>>>>>> ncl 328> rts at gsnYRefLine = 0. ;
>>>>>>>>>>> reference line
>>>>>>>>>>> ncl 329> rts at trXMaxF = max(lon)
>>>>>>>>>>> ncl 330> rts at trXMinF = min(lon)
>>>>>>>>>>> ncl 331>
>>>>>>>>>>> ncl 332> rts at pmLegendDisplayMode = "Always" ;
>>>>>>>>>>> turn on legend
>>>>>>>>>>> ncl 333> rts at pmLegendSide = "Top" ;
>>>>>>>>>>> Change location of
>>>>>>>>>>> ncl 334> rts at pmLegendParallelPosF = 1.16 ;
>>>>>>>>>>> move units right
>>>>>>>>>>> ncl 335> rts at pmLegendOrthogonalPosF = -0.50 ;
>>>>>>>>>>> move units down
>>>>>>>>>>> ncl 336> rts at pmLegendWidthF = 0.15 ;
>>>>>>>>>>> Change width and
>>>>>>>>>>> ncl 337> rts at pmLegendHeightF = 0.15 ;
>>>>>>>>>>> height of legend.
>>>>>>>>>>> ncl 338> rts at lgLabelFontHeightF = 0.0175
>>>>>>>>>>> ncl 339>
>>>>>>>>>>> ncl 340>
>>>>>>>>>>> ncl 341> rtsP = True ;
>>>>>>>>>>> modify the panel plot
>>>>>>>>>>> ncl 342> ; rtsP at gsnMaximize = True ;
>>>>>>>>>>> large format
>>>>>>>>>>> ncl 343> rtsP at gsnPanelMainString = "Multivariate EOF:
>>>>>>>>>>> 15S-15N: "+yrStrt+"-"+yrLast
>>>>>>>>>>> ncl 344>
>>>>>>>>>>> ncl 345> do n=0,neof-1
>>>>>>>>>>> ncl 346> rts at xyExplicitLegendLabels = (/"OLR:
>>>>>>>>>>> "+sprintf("%4.1f", pcv_eof_u200(n)) +"%" \
>>>>>>>>>>> ncl 346> ,"U850:
>>>>>>>>>>> "+sprintf("%4.1f", pcv_eof_u850(n))+"%" \
>>>>>>>>>>> ncl 346> ,"U200:
>>>>>>>>>>> "+sprintf("%4.1f", pcv_eof_olr(n))+"%" /)
>>>>>>>>>>> ncl 347> rts at gsnLeftString = "EOF "+(n+1)
>>>>>>>>>>> ncl 348> rts at gsnRightString = sprintf("%3.1f",ceof at pcvar(n))
>>>>>>>>>>> +"%"
>>>>>>>>>>> ncl 349> plot(n) = gsn_csm_xy (wks,lon,ceof(:,n,:),rts)
>>>>>>>>>>> ncl 350> end do
>>>>>>>>>>> ncl 351> gsn_panel(wks,plot(0:1),(/2,1/),rtsP) ; now draw
>>>>>>>>>>> as one plot
>>>>>>>>>>> ncl 352>
>>>>>>>>>>> ncl 353> ;-----------------------------------------
>>>>>>>>>>> ncl 354> ; The following doesn't work with some older versions
>>>>>>>>>>> of NCL
>>>>>>>>>>> ncl 355> ; With old versions, the user must delete each
>>>>>>>>>>> individually.
>>>>>>>>>>> ncl 356> ;-----------------------------------------
>>>>>>>>>>> ncl 357> delete([/ rts at xyExplicitLegendLabels,
>>>>>>>>>>> rts at pmLegendDisplayMode \
>>>>>>>>>>> ncl 357> , rts at xyLineThicknesses ,
>>>>>>>>>>> rts at gsnLeftString \
>>>>>>>>>>> ncl 357> , rts at gsnRightString ,
>>>>>>>>>>> rts at xyLineColors \
>>>>>>>>>>> ncl 357> , rts at trXMaxF , rts at trXMinF
>>>>>>>>>>> /] )
>>>>>>>>>>> ncl 358>
>>>>>>>>>>> ncl 359> lag = ispan(-mxlag,mxlag,1)
>>>>>>>>>>> ncl 360> lag at long_name = "lag (days)"
>>>>>>>>>>> ncl 361>
>>>>>>>>>>> ncl 362> plot(0) = gsn_csm_xy (wks, lag
>>>>>>>>>>> ,ccr_12,rts)
>>>>>>>>>>> ncl 363> rtsP at gsnPanelMainString = "Cross Correlation:
>>>>>>>>>>> Multivariate EOF: 15S-15N: " \
>>>>>>>>>>> ncl 363> + yrStrt+"-"+yrLast
>>>>>>>>>>> ncl 364> rtsP at gsnPaperOrientation = "portrait" ;
>>>>>>>>>>> force portrait
>>>>>>>>>>> ncl 365> gsn_panel(wks,plot(0),(/1,1/),rtsP) ; now draw as
>>>>>>>>>>> one plot
>>>>>>>>>>> ncl 366>
>>>>>>>>>>> ncl 367> ;************************************************
>>>>>>>>>>> ncl 368> ; MJO "strong" index
>>>>>>>>>>> ncl 369> ;************************************************
>>>>>>>>>>> ncl 370> rts at gsnYRefLine = 1.0
>>>>>>>>>>> ncl 371> rts at gsnYRefLineColor = "black"
>>>>>>>>>>> ncl 372> rts at xyMonoDashPattern = True
>>>>>>>>>>> ncl 373> rts at xyLineColors = (/"black", "blue"/)
>>>>>>>>>>> ncl 374> rts at xyLineThicknesses = (/1, 2/)
>>>>>>>>>>> ncl 375> rts at pmLegendDisplayMode = "Always" ;
>>>>>>>>>>> turn on legend
>>>>>>>>>>> ncl 376> rts at pmLegendWidthF = 0.12 ;
>>>>>>>>>>> Change width and
>>>>>>>>>>> ncl 377> rts at pmLegendHeightF = 0.10 ;
>>>>>>>>>>> height of legend.
>>>>>>>>>>> ncl 378> rts at pmLegendParallelPosF = 0.86 ;
>>>>>>>>>>> move units right
>>>>>>>>>>> ncl 379> rts at pmLegendOrthogonalPosF = -0.40 ;
>>>>>>>>>>> move units down
>>>>>>>>>>> ncl 380> rts at xyExplicitLegendLabels = (/"daily", "91-day
>>>>>>>>>>> runavg" /)
>>>>>>>>>>> ncl 381>
>>>>>>>>>>> ncl 382> mjo_ind_plt = new ( (/2,ntim/), typeof(mjo_ts_index))
>>>>>>>>>>> ncl 383> mjo_ind_plt(0,:) = mjo_ts_index
>>>>>>>>>>> ncl 384> mjo_ind_plt(1,:) = (/ mjo_ts_index_smt /)
>>>>>>>>>>> ncl 385> plot(0) = gsn_csm_xy(wks, yrfrac,mjo_ind_plt,rts)
>>>>>>>>>>> ncl 386>
>>>>>>>>>>> ncl 387> rtsP at gsnPanelMainString = "MJO Index: (PC1^2+
>>>>>>>>>>> PC2^2) : 15S-15N: "+yrStrt+"-"+yrLast
>>>>>>>>>>> ncl 388> gsn_panel(wks,plot(0),(/1,1/),rtsP) ; now draw as
>>>>>>>>>>> one plot
>>>>>>>>>>> ncl 389>
>>>>>>>>>>> ncl 390> end
>>>>>>>>>>>
>>>>>>>>>>> Thank you
>>>>>>>>>>>
>>>>>>>>>>> _______________________________________________
>>>>>>>>>>> ncl-talk mailing list
>>>>>>>>>>> ncl-talk at ucar.edu
>>>>>>>>>>> List instructions, subscriber options, unsubscribe:
>>>>>>>>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> --
>>>>>>>>>> Adam Phillips
>>>>>>>>>> Associate Scientist, Climate and Global Dynamics Laboratory,
>>>>>>>>>> NCAR
>>>>>>>>>> www.cgd.ucar.edu/staff/asphilli/ 303-497-1726
>>>>>>>>>>
>>>>>>>>>> <http://www.cgd.ucar.edu/staff/asphilli>
>>>>>>>>>>
>>>>>>>>> _______________________________________________
>>>>>>>>> ncl-talk mailing list
>>>>>>>>> ncl-talk at ucar.edu
>>>>>>>>> List instructions, subscriber options, unsubscribe:
>>>>>>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>>>>>>
>>>>>>>> _______________________________________________
>>>>> ncl-talk mailing list
>>>>> ncl-talk at ucar.edu
>>>>> List instructions, subscriber options, unsubscribe:
>>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>>
>>>>
>>>>
>>>> --
>>>> Barry H. Lynn, Ph.D
>>>> Senior Associate Scientist, Lecturer,
>>>> The Institute of the Earth Science,
>>>> The Hebrew University of Jerusalem,
>>>> Givat Ram, Jerusalem 91904, Israel
>>>> Tel: 972 547 231 170
>>>> Fax: (972)-25662581
>>>>
>>>> C.E.O, Weather It Is, LTD
>>>> Weather and Climate Focus
>>>> http://weather-it-is.com
>>>> Jerusalem, Israel
>>>> Local: 02 930 9525
>>>> Cell: 054 7 231 170
>>>> Int-IS: x972 2 930 9525
>>>>
>>>> _______________________________________________
>>> ncl-talk mailing list
>>> ncl-talk at ucar.edu
>>> List instructions, subscriber options, unsubscribe:
>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
>>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191223/58367aa4/attachment.html>
More information about the ncl-talk
mailing list