[ncl-talk] fatal:divide: Division by 0, Can't continuefatal:Div: operator failed, can't continue fatal:["Execute.c":8635]:Execute: Error occurred at or near line 129
Rahpeni Fajarianti
rahpenifajarianti at gmail.com
Wed Dec 18 08:13:41 MST 2019
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>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191218/a5243189/attachment.html>
More information about the ncl-talk
mailing list