[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
Fri Dec 20 12:21:40 MST 2019
[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@
> _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/20191220/3d445102/attachment.html>
More information about the ncl-talk
mailing list