[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

Barry Lynn barry.h.lynn at gmail.com
Thu Dec 19 23:15:28 MST 2019


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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191220/39fa97ee/attachment.html>


More information about the ncl-talk mailing list