<div dir="ltr">Hi Rahpeni,<div>Please always include ncl-talk on your replies; that way others can assist and see the correspondence. </div><div><br></div><div>You stated:</div><div><div>Thank you so much, it helps ! </div><div><br></div><div>I am still getting some errors.<br></div><div>I got :</div><div>fatal:Subscript out of range, error in subscript #1<br>fatal:An error occurred reading olr<br>fatal:["Execute.c":8635]:Execute: Error occurred at or near line 141</div></div><div><br></div><div><br></div><div>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):</div><div><br></div><div>ncl 140> Â  do ml=0,mlon-1<br>ncl 141> Â  Â  Â cdata(ml Â  Â  Â  ,:) = (/ Â olr(:,ml) /)<br>ncl 142> Â  Â  Â cdata(ml+ Â mlon,:) = (/ u850(:,ml) /)<br>ncl 143> Â  Â  Â cdata(ml+2*mlon,:) = (/ u200(:,ml) /)<br>ncl 144> Â  end do<br></div><div><br></div><div>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) /)</div><div><br></div><div>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. </div><div><br></div><div>Hope that helps!</div><div>Adam</div><div><br></div><div><br></div><div><br></div><div><br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Tue, Dec 17, 2019 at 1:40 AM Rahpeni Fajarianti via ncl-talk <<a href="mailto:ncl-talk@ucar.edu" target="_blank">ncl-talk@ucar.edu</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>Hi all ! <br></div><div>I want to make MJO clivar 14, but I have some problem.<br></div><div>This is :</div><div><div>fatal:divide: Division by 0, Can't continue</div>fatal:Div: operator failed, can't continue<br>fatal:["Execute.c":8635]:Execute: Error occurred at or near line 129</div><div><br></div><div>This is my script:<br></div><div>ncl 0> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/diagnostics_cam.ncl"</div>ncl 1> ;******************************************************<br> u850 Â = dim_rmvmean_n(u850, 0)<br>  Â u200 Â = dim_rmvmean_nncl 2> ;<br>(uncl 3> ; mjoclivar_14.ncl<br>ncl 4> ;<br>ncl 5> ;***********************************************************<br>ncl 6> ; Combined EOFs<br>hncl 7> ; Latest Update: July, 2016: Eun-Pa Lim; Bureau of Meteorology, Australia<br>ncl 8> ;***********************************************************<br>iancncl 9> ;;<br>nncl 10> ;; Â  Â  Â The following are automatically loaded from 6.2.0 onward<br>ncl 11> ;;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" Â <br>ncl 12> ;;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" Â <br>ncl 13> ;;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" <br>ncl 14> <br>ncl 15> undef("read_rename")<br>ncl 16> function read_rename(f[1]:file, varName[1]:string Â  Â  Â  \<br>ncl 16> Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  ,iStrt[1]:integer, iLast[1]:integer \<br>ncl 16> Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  ,latS[1]:numeric , latN[1]:numeric Â )<br>ncl 17> ; Utility to force specific named dimensions<br>ncl 18> ; This is done for historical reasons (convenience) <br>ncl 19> begin<br>ncl 20> Â  Â work Â  Â = f->$varName$(iStrt:iLast,{latS:latN},:) Â  ; (time,lat,lon)<br>ncl 21> Â  Â work!0 Â = "time" Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â ; CAM model names<br>ncl 22> Â  Â work!1 Â = "lat"<br>ncl 23> Â  Â work!2 Â = "lon"<br>ncl 24> Â  Â return(work)<br>ncl 25> end<br>ncl 26> ; =========================> Â MAIN Â <==============================<br>ncl 27> begin<br>ncl 28> Â  Â neof Â  Â = Â 2<br>ncl 29> <br>ncl 30> Â  Â latS Â  Â = -15<br>ncl 31> Â  Â latN Â  Â = Â 15<br>ncl 32> <br> ncl 33> Â  Â ymdStrt = 20180101 Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  ; start yyyymmdd<br>ncl 34> Â  Â ymdLast = 20181231 Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  ; last Â <br>ncl 35> <br>ncl 36> Â  Â yrStrt Â = ymdStrt/10000<br>ncl 37> Â  Â yrLast Â = ymdLast/10000<br>ncl 38> <br>ncl 39> Â  Â pltDir Â = "/home/peni/" Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  ; plot directory<br>*ncl 40> Â  Â pltType = "png" <br>ncl 41> Â  Â pltName = "mjoclivar" Â  <br>ncl 42> <br>,ncl 43> Â  Â diri Â  Â = "/home/peni/" Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  ; input directory<br>ncl 44> <br>ncl 45> Â  Â filolr Â = "<a href="http://olr.day.anomalies.2018.nc" target="_blank">olr.day.anomalies.2018.nc</a>"<br>ncl 46> Â  Â filu200 = "<a href="http://uwnd.anomalies.nc" target="_blank">uwnd.anomalies.nc</a>"<br>ncl 47> Â  Â filu850 = "<a href="http://uwnd2.anomalies.nc.nc" target="_blank">uwnd2.anomalies.nc.nc</a>"<br>ncl 48> <br>ncl 49> ;************************************************<br>ncl 50> ; create BandPass Filter<br>ncl 51> ;************************************************<br>ncl 52> Â  ihp Â  Â  Â = 2 Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  ; bpf=>band pass filter<br>ncl 53> Â  nWgt Â  Â  = 201<br>ncl 54> Â  sigma Â  Â = 1.0 Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  ; Lanczos sigma<br>ncl 55> Â  fca Â  Â  Â = 1./100.<br>ncl 56> Â  fcb Â  Â  Â = 1./20.<br>ncl 57> Â  wgt Â  Â  Â = filwgts_lanczos (nWgt, ihp, fca, fcb, sigma )<br>ncl 58> <br>)ncl 59> ;***********************************************************<br>ncl 60> ; Find the indices corresponding to the start/end times<br>ncl 61> ;***********************************************************<br>ncl 62> Â  Â f Â  Â  Â  = addfile (diri+filolr , "r") Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  <br>ncl 63> Â  Â TIME Â  Â = f->time Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â ; days since ...<br>ncl 64> <br> ncl 65> Â  Â YMD Â  Â  = cd_calendar(TIME, -2) Â  Â  Â  Â  Â  Â ; entire (time,6)<br>ncl 66> <br>ncl 67> Â  Â iStrt Â  = ind(YMD.eq.ymdStrt) Â  Â  Â  Â  Â  Â  Â ; index start<br>ncl 68> Â  Â iLast Â  = ind(YMD.eq.ymdLast) Â  Â  Â  Â  Â  Â  Â ; index last <br>ncl 69> Â  Â delete([/ TIME, YMD /])<br>ncl 70> <br>ncl 71> ;***********************************************************<br>ncl 72> ; Read anomalies<br>ncl 73> ;***********************************************************<br>ncl 74> <br>ncl 75> Â  Â work Â  Â = read_rename(f,"OLR_anom",iStrt,iLast,latS,latN) ; (time,lat,lon) <br>ncl 76> Â  Â OLR Â  Â  = dim_avg_n_Wrap(work, 1) Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  ; (time,lon)<br>ncl 77> Â  Â delete(work)<br>ncl 78> Â  Â <br>ncl 79> Â  Â f Â  Â  Â  = addfile (diri+filu850 , "r") Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  <br>ncl 80> Â  Â work Â  Â = read_rename(f,"U_anom",iStrt,iLast,latS,latN) ; (time,lat,lon) <br>ncl 81> Â  Â U850 Â  Â = dim_avg_n_Wrap(work, 1) Â  Â  Â  Â  Â ; (time,lon)<br>ncl 82> Â  Â delete(work)<br>ncl 83> Â  Â <br>ncl 84> Â  Â f Â  Â  Â  = addfile (diri+filu200 , "r") Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  <br>ncl 85> Â  Â work Â  Â = read_rename(f,"U_anom",iStrt,iLast,latS,latN) ; (time,lat,lon) <br>ncl 86> Â  Â U200 Â  Â = dim_avg_n_Wrap(work, 1) Â  Â  Â  Â  Â ; (time,lon)<br>ncl 87> <br>ncl 88> Â  Â dimw Â  Â = dimsizes( work )<br>ncl 89> Â  Â ntim Â  Â = dimw(0)<br>ncl 90> Â  Â nlat Â  Â = dimw(1)<br>ncl 91> Â  Â mlon Â  Â = dimw(2)<br>ncl 92> Â  Â delete(work)<br>ncl 93> <br>ncl 94> Â  Â lon Â  Â  = OLR&lon Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â <br>ncl 95> Â  Â time Â  Â = OLR&time Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  <br>ncl 96> Â  Â date Â  Â = cd_calendar(time, -2) Â  Â  Â  Â  Â  Â ; yyyymmdd<br>ncl 97> <br>ncl 98> ;************************************************<br>ncl 99> ; Apply the band pass filter to the original anomalies<br>ncl 100> ;************************************************<br>ncl 101> Â  Â olr Â  = wgt_runave_n_Wrap ( OLR, wgt, 0, 0) ; (time,lon)<br>) <br>;************************************************<br>  Â  Â  Â  Â  Â  Â  Â  Â  Â <br>  imax_olr_eof1 Â  = maxind(ceof(0,0,:)) Â  Â  Â  Â  <br>  imax_olr_eof2 ncl 102> Â  Â u850 Â = wgt_runave_n_Wrap (U850, wgt, 0, 0)<br>ncl 103> Â  Â u200 Â = wgt_runave_n_Wrap (U200, wgt, 0, 0)<br>ncl 104> <br>ncl 105> ;************************************************<br>ncl 106> ; remove temporal means of band pass series: *not* necessary <br>ncl 107> ;************************************************<br>ncl 108> Â  Â olr Â  = dim_rmvmean_n( olr, 0) Â  Â  Â  Â  Â  Â  Â ; (time,lon)<br>ncl 109> Â  Â u850 Â = dim_rmvmean_n(u850, 0)<br>ncl 110> Â  Â u200 Â = dim_rmvmean_n(u200, 0)<br>ncl 111> <br>ncl 112> ;************************************************<br>ncl 113> ; Compute the temporal variance at each lon<br>ncl 114> ;************************************************<br>ncl 115> Â  Â var_olr Â = dim_variance_n_Wrap( olr, 0) Â  Â  ; (lon)<br>ncl 116> Â  Â var_u850 = dim_variance_n_Wrap(u850, 0)<br>ncl 117> Â  Â var_u200 = dim_variance_n_Wrap(u200, 0)<br>ncl 118> <br>ncl 119> ;************************************************<br>ncl 120> ; Compute the zonal mean of the temporal variance<br> Oncl 121> ;************************************************<br>ncl 122> Â  zavg_var_olr Â = dim_avg_n_Wrap( var_olr , 0) Â  Â  <br>ncl 123> Â  zavg_var_u850 = dim_avg_n_Wrap( var_u850, 0)<br>ncl 124> Â  zavg_var_u200 = dim_avg_n_Wrap( var_u200, 0)<br>ncl 125> <br>ncl 126> ;************************************************<br>ncl 127> ; Normalize by sqrt(avg_var*)<br>ncl 128> ;************************************************<br>ncl 129> Â  olr Â  = Â olr/sqrt(zavg_var_olr ) Â  Â  Â  Â  Â ; (time,lon)<br>ncl 130> Â  u850 Â = u850/sqrt(zavg_var_u850)<br>ncl 131> Â  u200 Â = u200/sqrt(zavg_var_u200)<br>ncl 132> <br>cncl 133> ;************************************************<br>ncl 134> ; Combine the normalized data into one variable<br>ncl 135> ;************************************************<br>ncl 136> Â  cdata Â  Â  = new ( (/3*mlon,ntim/), typeof(olr), getFillValue(olr))<br>ncl 137> Â  do ml=0,mlon-1<br>ncl 138> Â  Â  Â cdata(ml Â  Â  Â  ,:) = (/ Â olr(:,ml) /)<br>ncl 139> Â  Â  Â cdata(ml+ Â mlon,:) = (/ u850(:,ml) /)<br>ncl 140> Â  Â  Â cdata(ml+2*mlon,:) = (/ u200(:,ml) /)<br>ncl 141> Â  end do<br>ncl 142> <br>=ncl 143> ;************************************************<br>ncl 144> ; Compute **combined** EOF; Sign of EOF is arbitrary<br>ncl 145> ;************************************************<br>ncl 146> Â  eof_cdata Â  Â = eofunc(cdata Â  , neof, False) Â  Â  Â ; (neof,3*mlon)<br>ncl 147> Â  print("==============")<br>ncl 148> Â  printVarSummary(eof_cdata)<br>ncl 149> Â  printMinMax(eof_cdata, True)<br>ncl 150> <br>ncl 151> Â  eof_ts_cdata = eofunc_ts(cdata,eof_cdata,False) Â  ; (neof,3*ntim)<br>ncl 152> Â  print("==============") Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â <br>ncl 153> Â  printVarSummary(eof_ts_cdata)<br>ncl 154> Â  printMinMax(eof_ts_cdata, True)<br>ncl 155> <br>incl 156> ;************************************************<br>ncl 157> ; For clarity, explicitly extract each variable. Create time series <br>ncl 158> ;************************************************<br>ncl 159> <br>tncl 160> Â  nvar = 3 Â ; "olr", "u850", "u200"<br>ncl 161> Â  ceof = new( (/nvar,neof,mlon/), typeof(cdata), getFillValue(cdata))<br>ncl 162> <br>dncl 163> Â  do n=0,neof-1<br>ncl 164> Â  Â  Â ceof(0,n,:) = eof_cdata(n,0:mlon-1) Â  Â  Â ; olr<br>ncl 165> Â  Â  Â ceof(1,n,:) = eof_cdata(n,mlon:2*mlon-1) ; u850<br>ncl 166> Â  Â  Â ceof(2,n,:) = eof_cdata(n,2*mlon:) Â  Â  Â  ; u200<br>ncl 167> Â  end do<br>ncl 168> <br>(ncl 169> Â  ceof!0 Â  = "var"<br>ncl 170> Â  ceof!1 Â  = "eof"<br>ncl 171> Â  ceof!2 Â  = "lon" Â  <br>ncl 172> Â  ceof&lon = Â olr&lon<br>ncl 173> <br>incl 174> Â  ceof_ts Â  Â  Â  Â = new( (/nvar,neof,ntim/), typeof(cdata), getFillValue(cdata))<br>ncl 175> Â  ceof_ts(0,:,:) = eofunc_ts_Wrap( olr(lon|:,time|:),ceof(0,:,:),False) Â  ; (0,neof,ntim)<br>ncl 176> Â  ceof_ts(1,:,:) = eofunc_ts_Wrap(u850(lon|:,time|:),ceof(1,:,:),False) Â  ; (1,neof,ntim)<br>ncl 177> Â  ceof_ts(2,:,:) = eofunc_ts_Wrap(u200(lon|:,time|:),ceof(2,:,:),False) Â  ; (2,neof,ntim)<br>ncl 178> <br>ncl 179> ;**********************************************t*<br>ncl 180> ; Add code contributed by Marcus N. Morgan, Florida Institute of Technology; Feb 2015<br>ncl 181> ; Calculate % variance (pcv_ )accounted for by OLR, U850 and U200<br>ncl 182> ;************************************************<br>ncl 183> <br>ncl 184> Â  Â  pcv_eof_olr Â = new(neof,typeof(ceof))<br>ncl 185> Â  Â  pcv_eof_u850 = new(neof,typeof(ceof))<br>ncl 186> Â  Â  pcv_eof_u200 = new(neof,typeof(ceof))<br>ncl 187> Â  Â  Â  <br>ncl 188> Â  Â  do n=0,neof-1<br>ncl 189> Â  Â  Â  Â pcv_eof_olr(n) Â = avg((ceof(0,n,:)*sqrt(ceof@eval(n)))^2)*100<br>ncl 190> Â  Â  Â  Â pcv_eof_u850(n) = avg((ceof(1,n,:)*sqrt(ceof@eval(n)))^2)*100<br>ncl 191> Â  Â  Â  Â pcv_eof_u200(n) = avg((ceof(2,n,:)*sqrt(ceof@eval(n)))^2)*100<br>ncl 192> Â  Â  Â ;;print("pcv: neof="+(n+1)+": Â "+pcv_eof_olr(n)+" Â "+pcv_eof_u850(n)+" Â "+pcv_eof_u200(n))<br>ncl 193> Â  Â  end do<br>ncl 194> <br> ncl 195> ;************************************************<br>ncl 196> ; Change sign of EOFs for spatial structures of PC1 and PC2 <br>ncl 197> ; to represent convection over the tropical Indian Ocean and the tropical western Pacific Ocean, respectively <br>ncl 198> ; (Ad hoc approach) <br>ncl 199> ;************************************************<br>ncl 200> Â  Â  Â  Â  Â  Â  Â  Â  Â  Â <br>ncl 201> Â  imax_olr_eof1 Â  = maxind(ceof(0,0,:)) Â  Â  Â  Â  <br>ncl 202> Â  imax_olr_eof2 Â  = maxind(ceof(0,1,:)) Â  Â <br>ncl 203> <br>)ncl 204> Â  lonmax_eof1 = ceof&lon(imax_olr_eof1) Â  Â  Â ; longitude of max value (i.e. suppressed convection)<br>ncl 205> Â  lonmax_eof2 = ceof&lon(imax_olr_eof2)<br>ncl 206> <br>sncl 207> Â  if (lonmax_eof1.ge.60 .and. lonmax_eof1.lt.180) then Â ; Change the sign of EOF1 <br>ncl 208> Â  Â  Â  ceof(:,0,:) Â  Â  Â  = -ceof(:,0,:) Â  Â  Â  Â  Â  Â  Â  Â  Â ; if OLR is positive<br>ncl 209> Â  Â  Â  ceof_ts(:,0,:) Â  Â = -ceof_ts(:,0,:) Â  Â  Â  Â  Â  Â  Â  ; Â over the tropical Indian Ocean<br>ncl 210> Â  Â  Â  eof_cdata(0,:) Â  Â = -eof_cdata(0,:)<br>ncl 211> Â  Â  Â  eof_ts_cdata(0,:) = -eof_ts_cdata(0,:)<br>ncl 212> Â  end if<br>ncl 213> <br>encl 214> Â  if (lonmax_eof2.ge.120 .and. lonmax_eof2.lt.180) then Â ; Change the sign of EOF2<br>ncl 215> Â  Â  Â  ceof(:,1,:) Â  Â  Â  = -ceof(:,1,:) Â  Â  Â  Â  Â  Â  Â  Â  Â  ; if OLR is positive <br>ncl 216> Â  Â  Â  ceof_ts(:,1,:) Â  Â = -ceof_ts(:,1,:) Â  Â  Â  Â  Â  Â  Â  Â ; over the tropical western Pacific Ocean<br>ncl 217> Â  Â  Â  eof_cdata(1,:) Â  Â = -eof_cdata(1,:)<br>ncl 218> Â  Â  Â  eof_ts_cdata(1,:) = -eof_ts_cdata(1,:)<br>ncl 219> Â  end if<br>ncl 220> <br>ncl 221> Â  print("==============")<br>rncl 222> Â  printVarSummary(eof_cdata)<br>ncl 223> Â  printMinMax(eof_cdata, True)<br>ncl 224> <br>ncl 225> ;************************************************<br>ncl 226> ; Compute cross correlation of each variable's EOF time series at zero-lag<br>ncl 227> ;************************************************<br>ncl 228> Â  r_olr_u850 Â = escorc(ceof_ts(0,:,:) , ceof_ts(1,:,:) ) Â ; (neof)<br>ncl 229> Â  r_olr_u200 Â = escorc(ceof_ts(0,:,:) , ceof_ts(2,:,:) )<br>ncl 230> Â  r_u850_u200 = escorc(ceof_ts(1,:,:) , ceof_ts(2,:,:) )<br>ncl 231> <br>sncl 232> Â  print("==============")<br>ncl 233> Â  do n=0,neof-1<br>  ncl 234> Â  Â  Â print("neof="+n \<br>ncl 234> Â  Â  Â  Â  Â  +" Â r_olr_u850=" +sprintf("%4.3f",r_olr_u850(n)) Â \<br>ncl 234> Â  Â  Â  Â  Â  +" Â r_olr_u200=" +sprintf("%4.3f",r_olr_u200(n)) Â \<br>ncl 234> Â  Â  Â  Â  Â  +" Â r_u850_u200="+sprintf("%4.3f",r_u850_u200(n)) )<br>ncl 235> Â  end do<br>ncl 236> Â  print("==============")<br>ncl 237> <br>ncl 238> ;************************************************<br>5N: "+yrStrt+"-"+yrLast <br>  <br>  do ncl 239> ; Compute cross correlation of the multivariate EOF; EOF 1 vs EOF 2<br>ncl 240> ;************************************************<br>ncl 241> <br> ncl 242> Â  mxlag Â  Â  = 25<br>ncl 243> Â  rlag_01 Â  = esccr(eof_ts_cdata(0,:),eof_ts_cdata(1,:), mxlag) Â  ; (N,mxlag+1)<br>ncl 244> Â  rlag_10 Â  = esccr(eof_ts_cdata(1,:),eof_ts_cdata(0,:), mxlag) Â  ; (N,mxlag+1)<br>ncl 245> Â  ccr_12 Â  Â = new ( (/2*mxlag+1/), float) Â  Â <br>ncl 246> <br>cncl 247> Â  ccr_12(mxlag:) Â  Â = rlag_10(0:mxlag) Â  <br>ncl 248> Â  ccr_12(0:mxlag) Â  = rlag_01(::-1) Â  Â  Â  ; reverse order<br>ncl 249> ;;print(ccr_12)<br>ncl 250> <br>sncl 251> <br>Pncl 252> ;************************************************<br>ncl 253> ; Normalize the multivariate EOF 1&2 component time series<br>ncl 254> ; Compute (PC1^2+PC2^2): values > 1 indicate "strong" periods<br>ncl 255> ;************************************************<br>ncl 256> Â  eof_ts_cdata(0,:) = eof_ts_cdata(0,:)/stddev(eof_ts_cdata(0,:))<br>ncl 257> Â  eof_ts_cdata(1,:) = eof_ts_cdata(1,:)/stddev(eof_ts_cdata(1,:))<br>ncl 258> <br>ncl 259> Â  mjo_ts_index Â  Â  Â = eof_ts_cdata(0,:)^2 + eof_ts_cdata(1,:)^2 <br>ncl 260> Â  mjo_ts_index_smt Â = runave(mjo_ts_index, 91, 0) ; 91-day running mean<br>ncl 261> <br> ncl 262> Â  nGood Â  = num(.not.ismissing(mjo_ts_index)) Â  Â  ; # non-missing<br>ncl 263> Â  nStrong = num(mjo_ts_index .ge. 1.0)<br>ncl 264> Â  print("nGood="+nGood+" Â  nStrong="+nStrong+" Â  nOther="+(nGood-nStrong))<br>ncl 265> <br>rncl 266> ;************************************************<br>ncl 267> ; Write PC results to netCDF for use in another example.<br>ncl 268> ;************************************************<br>ncl 269> Â  mjo_ts_index!0 Â  Â = "time"<br>ncl 270> Â  mjo_ts_index&time = time <br>ncl 271> Â  mjo_ts_index@long_name = "MJO PC INDEX" <br>ncl 272> Â  mjo_ts_index@info Â  Â  Â = "(PC1^2 + PC2^2)" <br>ncl 273> <br> ncl 274> Â  PC1 Â  Â  Â  Â  Â  = eof_ts_cdata(0,:)<br>ncl 275> Â  PC1!0 Â  Â  Â  Â  = "time"<br>ncl 276> Â  PC1&time Â  Â  Â = Â time<br>ncl 277> Â  PC1@long_name = "PC1"<br>ncl 278> Â  PC1@info Â  Â  Â = "PC1/stddev(PC1)"<br>ncl 279> <br>oncl 280> Â  PC2 Â  Â  Â  Â  Â  = eof_ts_cdata(1,:)<br>ncl 281> Â  PC2!0 Â  Â  Â  Â  = "time"<br>ncl 282> Â  PC2&time Â  Â  Â = Â time<br>ncl 283> Â  PC2@long_name = "PC2"<br>ncl 284> Â  PC2@info Â  Â  Â = "PC2/stddev(PC2)"<br>ncl 285> <br>ncl 286> Â  diro = "./"<br>ncl 287> Â  filo = "MJO_PC_INDEX.nc"<br>ncl 288> Â  system("/bin/rm -f "+diro+filo) Â  ; remove any pre-existing file<br>ncl 289> Â  ncdf = addfile(diro+filo,"c") Â  Â  ; open output netCDF file<br>of legend.<br>  ncl 290> Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  ; make time an UNLIMITED dimension <br>ncl 291> Â  filedimdef(ncdf,"time",-1,True) Â  ; recommended Â for most applications<br>ncl 292> Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  ; output variables directly<br>ncl 293> Â  ncdf->MJO_INDEX = mjo_ts_index Â  Â <br>ncl 294> Â  ncdf->PC1 Â  Â  Â  = PC1 Â  Â  <br>ncl 295> Â  ncdf->PC2 Â  Â  Â  = PC2 Â  Â  <br>ncl 296> <br>ncl 297> ;------------------------------------------------------------<br>ncl 298> ; PLOTS<br>ncl 299> ;------------------------------------------------------------<br>ncl 300> <br>ncl 301> Â  yyyymmdd = cd_calendar(time, -2)<br>tncl 302> Â  yrfrac Â  = yyyymmdd_to_yyyyfrac(yyyymmdd, 0.0)<br>ncl 303> Â  delete([/ yrfrac@long_name, lon@long_name /])<br>ncl 304> <br>ncl 305> Â  day Â  Â  Â = ispan(-mxlag, mxlag, 1)<br>ncl 306> Â ;day@long_name = "lag (day)"<br>ncl 307> <br>ncl 308> Â  pltPath = pltDir+pltName<br>ncl 309> <br>ncl 310> Â  wks = gsn_open_wks(pltType,pltPath)<br>ncl 311> Â  gsn_define_colormap(wks,"default")<br>ncl 312> Â  plot = new(3,graphic) Â  Â  Â  Â  Â  Â  Â  Â <br>ncl 313> <br>ncl 314> ;************************************************<br>ncl 315> ; Multivariate EOF plots<br>ncl 316> ;************************************************<br>ncl 317> Â  rts Â  Â  Â  Â  Â  = True<br>ncl 318> Â  rts@gsnDraw Â  = False Â  Â  Â  ; don't draw yet<br>ncl 319> Â  rts@gsnFrame Â = False Â  Â  Â  ; don't advance frame yet<br>ncl 320> Â  rts@gsnScale Â = True Â  Â  Â  Â ; force text scaling Â  Â  Â  Â  Â  Â  Â  <br>ncl 321> <br>ncl 322> Â  rts@vpHeightF = 0.40 Â  Â  Â  Â ; Changes the aspect ratio<br>ncl 323> Â  rts@vpWidthF Â = 0.85<br>ncl 324> Â  rts@vpXF Â  Â  Â = 0.10 Â  Â  Â  Â ; change start locations<br>ncl 325> Â  rts@vpYF Â  Â  Â = 0.75 Â  Â  Â  Â ; the plot<br>ncl 326> Â  rts@xyLineThicknesses = (/2, 2, 2/)<br>ncl 327> Â  rts@xyLineColors Â  Â  Â = (/"black","red","blue"/)<br>ncl 328> Â  rts@gsnYRefLine Â  Â  Â  = 0. Â  Â  Â  Â  Â  Â  Â  Â  Â ; reference line Â  <br>ncl 329> Â  rts@trXMaxF Â  Â  Â  Â  Â  = max(lon)<br>ncl 330> Â  rts@trXMinF Â  Â  Â  Â  Â  = min(lon)<br>ncl 331> <br>ncl 332> Â  rts@pmLegendDisplayMode Â  Â = "Always" Â  Â  Â  Â  Â  Â ; turn on legend<br>ncl 333> Â  rts@pmLegendSide Â  Â  Â  Â  Â  = "Top" Â  Â  Â  Â  Â  Â  Â  ; Change location of <br>ncl 334> Â  rts@pmLegendParallelPosF Â  = 1.16 Â  Â  Â  Â  Â  Â  Â  Â ; move units right<br>ncl 335> Â  rts@pmLegendOrthogonalPosF = -0.50 Â  Â  Â  Â  Â  Â  Â  ; move units down<br>ncl 336> Â  rts@pmLegendWidthF Â  Â  Â  Â  = 0.15 Â  Â  Â  Â  Â  Â  Â  Â ; Change width and<br>ncl 337> Â  rts@pmLegendHeightF Â  Â  Â  Â = 0.15 Â  Â  Â  Â  Â  Â  Â  Â ; height of legend.<br>ncl 338> Â  rts@lgLabelFontHeightF Â  Â  = 0.0175<br>ncl 339> <br>ncl 340> <br>ncl 341> Â  rtsP Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  = True Â  Â  Â  Â  Â  Â  Â  Â ; modify the panel plot<br>ncl 342> ; Â rtsP@gsnMaximize Â  Â  Â  Â  Â  = True Â  Â  Â  Â  Â  Â  Â  Â ; large format<br>ncl 343> Â  rtsP@gsnPanelMainString Â  Â  = "Multivariate EOF: 15S-15N: "+yrStrt+"-"+yrLast <br>ncl 344> Â  <br>ncl 345> Â  do n=0,neof-1<br>ncl 346> Â  Â  rts@xyExplicitLegendLabels = (/"OLR: "+sprintf("%4.1f", pcv_eof_u200(n)) +"%" \<br>ncl 346> Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  ,"U850: "+sprintf("%4.1f", pcv_eof_u850(n))+"%" \<br>ncl 346> Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  ,"U200: "+sprintf("%4.1f", pcv_eof_olr(n))+"%" /)<br>ncl 347> Â  Â  rts@gsnLeftString Â = "EOF "+(n+1)<br>ncl 348> Â  Â  rts@gsnRightString = sprintf("%3.1f",ceof@pcvar(n)) Â +"%"<br>ncl 349> Â  Â  plot(n) = gsn_csm_xy (wks,lon,ceof(:,n,:),rts)<br>ncl 350> Â  end do<br>ncl 351> Â  gsn_panel(wks,plot(0:1),(/2,1/),rtsP) Â  Â  ; now draw as one plot<br>ncl 352> <br>ncl 353> ;-----------------------------------------<br>ncl 354> ; The following doesn't work with some older versions of NCL<br>ncl 355> ; With old versions, the user must delete each individually.<br>ncl 356> ;----------------------------------------- <br>ncl 357> Â  delete([/ rts@xyExplicitLegendLabels, rts@pmLegendDisplayMode \<br>ncl 357> Â  Â  Â  Â  Â  , rts@xyLineThicknesses Â  Â  , rts@gsnLeftString Â  Â  Â  \<br>ncl 357> Â  Â  Â  Â  Â  , rts@gsnRightString Â  Â  Â  Â , rts@xyLineColors Â  Â  Â  Â \<br>ncl 357> Â  Â  Â  Â  Â  , rts@trXMaxF Â  Â  Â  Â  Â  Â  Â  , rts@trXMinF Â  Â  Â  Â  Â  Â  /] )<br>ncl 358> <br>ncl 359> Â  lag Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â = ispan(-mxlag,mxlag,1)<br>ncl 360> Â  lag@long_name Â  Â  Â  Â  Â  Â  Â = "lag (days)"<br>ncl 361> <br>ncl 362> Â  plot(0) Â  Â  Â  Â  Â  Â  Â  Â  Â  Â = gsn_csm_xy (wks, lag ,ccr_12,rts)<br>ncl 363> Â  rtsP@gsnPanelMainString Â  Â = "Cross Correlation: Multivariate EOF: 15S-15N: " \<br>ncl 363> Â  Â  Â  Â  Â  Â  Â  Â  Â  Â + Â yrStrt+"-"+yrLast <br>ncl 364> Â  rtsP@gsnPaperOrientation Â  = "portrait" Â  Â  Â  Â ; force portrait<br>ncl 365> Â  gsn_panel(wks,plot(0),(/1,1/),rtsP) Â  Â  ; now draw as one plot<br>ncl 366> <br>ncl 367> ;************************************************<br>ncl 368> ; MJO "strong" index <br>ncl 369> ;************************************************<br>ncl 370> Â  rts@gsnYRefLine Â  Â  Â  Â = 1.0<br>ncl 371> Â  rts@gsnYRefLineColor Â  = "black"<br>ncl 372> Â  rts@xyMonoDashPattern Â = True<br>ncl 373> Â  rts@xyLineColors Â  Â  Â  = (/"black", "blue"/)<br>ncl 374> Â  rts@xyLineThicknesses Â = (/1, 2/)<br>ncl 375> Â  rts@pmLegendDisplayMode Â  Â = "Always" Â  Â  Â  Â  Â  Â ; turn on legend<br>ncl 376> Â  rts@pmLegendWidthF Â  Â  Â  Â  = 0.12 Â  Â  Â  Â  Â  Â  Â  Â ; Change width and<br>ncl 377> Â  rts@pmLegendHeightF Â  Â  Â  Â = 0.10 Â  Â  Â  Â  Â  Â  Â  Â ; height of legend.<br>ncl 378> Â  rts@pmLegendParallelPosF Â  = 0.86 Â  Â  Â  Â  Â  Â  Â  Â ; move units right<br>ncl 379> Â  rts@pmLegendOrthogonalPosF = -0.40 Â  Â  Â  Â  Â  Â  Â  ; move units down<br>ncl 380> Â  rts@xyExplicitLegendLabels = (/"daily", "91-day runavg" /)<br>ncl 381> <br>ncl 382> Â  mjo_ind_plt = new ( (/2,ntim/), typeof(mjo_ts_index))<br>ncl 383> Â  mjo_ind_plt(0,:) = mjo_ts_index<br>ncl 384> Â  mjo_ind_plt(1,:) = (/ mjo_ts_index_smt /)<br>ncl 385> Â  plot(0) = gsn_csm_xy(wks, yrfrac,mjo_ind_plt,rts)<br>ncl 386> <br>ncl 387> Â  rtsP@gsnPanelMainString Â  = "MJO Index: (PC1^2+ PC2^2) : 15S-15N: "+yrStrt+"-"+yrLast <br>ncl 388> Â  gsn_panel(wks,plot(0),(/1,1/),rtsP) Â  Â  ; now draw as one plot<br>ncl 389> <br>ncl 390> Â end<br><div><br></div><div>Thank you<br></div><br></div>
_______________________________________________<br>
ncl-talk mailing list<br>
<a href="mailto:ncl-talk@ucar.edu" target="_blank">ncl-talk@ucar.edu</a><br>
List instructions, subscriber options, unsubscribe:<br>
<a href="http://mailman.ucar.edu/mailman/listinfo/ncl-talk" rel="noreferrer" target="_blank">http://mailman.ucar.edu/mailman/listinfo/ncl-talk</a></blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div><div><div><span><font color="#888888">Adam Phillips <br></font></span></div><span><font color="#888888">Associate Scientist,  </font></span><span><font color="#888888">Climate and Global Dynamics Laboratory, NCAR<br></font></span></div></div><div><span><font color="#888888"><a href="http://www.cgd.ucar.edu/staff/asphilli/" target="_blank">www.cgd.ucar.edu/staff/asphilli/</a>   </font></span><span><font color="#888888">303-497-1726 </font></span></div><span><font color="#888888"></font></span><div><div><span><font color="#888888"><br></font></span><div><span><font color="#888888"><a href="http://www.cgd.ucar.edu/staff/asphilli" target="_blank"></a></font></span></div></div></div></div></div></div></div></div></div></div></div>