[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
Rahpeni Fajarianti
rahpenifajarianti at gmail.com
Fri Dec 20 08:12:57 MST 2019
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@
_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
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191220/71c1b05e/attachment-0001.html>
More information about the ncl-talk
mailing list