[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 16:57:38 MST 2019
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
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191219/ce8b620b/attachment.html>
More information about the ncl-talk
mailing list