[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
Sun Dec 29 01:45:58 MST 2019


Thank you for your advice. I am trying on it now

On Mon, Dec 23, 2019 at 10:32 PM Dennis Shea <shea at ucar.edu> wrote:

> No!
> ---
> mjoclivar_14 *requires multiple years of daily anomaly data*. Note the
> example uses anomalies from:
>
>    ymdStrt = 19800101                         ; start yyyymmdd
>    ymdLast = 20051231                         ; last
>
> Note the Lanczos band-pass filter specified by the MJO Clivar diagnostics. See top of section.
> ;************************************************
> ; create BandPass Filter
> ;************************************************
>   ihp      = 2                             ; bpf=>band pass filter
>   *nWgt     = 201*
>   sigma    = 1.0                           ; Lanczos sigma
>   fca      = 1./100.
>   fcb      = 1./20.
>   wgt      = *filwgts_lanczos* <http://www.ncl.ucar.edu/Document/Functions/Built-in/filwgts_lanczos.shtml> (nWgt, ihp, fca, fcb, sigma )
>
> *Hence to capture year 2018 you need (nWgt/2=201/2=) 150 or more days into year 2019.*
>
> ---
>
> *After* doing mjoclivar_14 for multiple years, you then specifya period in mjoclivar_15.
>
>
> On Fri, Dec 20, 2019 at 7:23 PM Rahpeni Fajarianti <
> rahpenifajarianti at gmail.com> wrote:
>
>> Thank you.
>>
>> Now, I use long term means for daily values for 1981-2010 as base
>> climatology data. It works.
>> I got min max values for :
>> - OLR anomaly (min=-185.5106410148043   max=77.65257890949087)
>> - 200hPa zonal wind anomaly  (min=-48.11269687686956
>> max=57.21568789939967)
>> - 850hPa zonal wind anomaly (min=-31.22578327383305
>> max=41.24650545117836).
>>
>> then, I use mjoclivar_14.ncl to get multivariate eof for data 2018. I
>> attach my script.
>>
>> Some error shown like this:
>> 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 converge
>> warning:eofunc_ts: 2 eigenvectors failed to converge
>> warning:eofunc_ts: 2 eigenvectors failed to converge
>> warning:avg: Entire input array contains missing values, can't compute
>> average
>>
>> fatal:Subscript out of range, error in subscript #0
>> fatal:An error occurred reading unknown
>> fatal:["Execute.c":8635]:Execute: Error occurred at or near line 207
>>
>> this is the error line :
>> 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)
>>
>> The data I used:
>>  olranom.nc
>> <https://drive.google.com/file/d/1xbmRiRO9VrTQ6qqKTfTq_rraiRSZFOXg/view?usp=drive_web>
>>  uanom.200.nc
>> <https://drive.google.com/file/d/1u1IYOGmjLuWAR4izwLkAgEMzWDxrkZlX/view?usp=drive_web>
>>  uanom.850.nc
>> <https://drive.google.com/file/d/1ZTtEukGralmUuwBWLGJ9ei6-N0yf199b/view?usp=drive_web>
>> The script I used:
>>
>> On Sat, Dec 21, 2019 at 2:21 AM Dennis Shea <shea at ucar.edu> wrote:
>>
>>>  [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 at _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/20191229/c413639b/attachment.html>


More information about the ncl-talk mailing list