<div dir="ltr"><div> [1] <br></div><div>PLEASE do not include your  line-by-line script in you emails. It makes the emails too long. Attach the script.</div><div><br></div><div>[2]</div><div>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.<br></div><div><br></div><div>[3]</div><div>I speculate you are using daily values from 2018 only. That is ***NOT*** a climatology. [Look up definition of climatology.] <br></div><div>There are no anomalies when you use one year of data That is why min-anomaly=0.0 md max-anomaly=0.0 </div><div><br></div><div><br></div><div><br></div><div><br></div><div> <br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Fri, Dec 20, 2019 at 10:14 AM Rahpeni Fajarianti via ncl-talk <<a href="mailto:ncl-talk@ucar.edu" target="_blank">ncl-talk@ucar.edu</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>I did it. Thanks. <br></div><div><br></div><div>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. <br></div><div><br></div><div><b>This is the result when I use printMinMax statement :</b> <br></div><div><br></div><div>Variable: olr<br>Type: double<br>Total Size: 5504200 bytes<br>            688025 values<br>Number of Dimensions: 3<br>Dimensions and sizes:  [time | 365] x [lat | 13] x [lon | 145]<br>Coordinates: <br>            time: [   0..524160]<br>            lat: [ -15..  15]<br>            lon: [   0.. 360]<br>Number Of Attributes: 1<br>  _FillValue :   -999000000<br>ncl 19>   printMinMax(olr,0)<br><b>(0)      min=-185.5106410148043   max=77.65257890949087</b><br>ncl 20>   print("------------------------------")<br>(0) ------------------------------<br>ncl 21>   printVarSummary(u200)<br><br>Variable: u200<br>Type: double<br>Total Size: 5504200 bytes<br>            688025 values<br>Number of Dimensions: 3<br>Dimensions and sizes: [time | 365] x [lat | 13] x [lon | 145]<br>Coordinates: <br>            time: [   0..524160]<br>            lat: [ -15..  15]<br>            lon: [   0.. 360]<br>Number Of Attributes: 1<br>  _FillValue :   -999000000<br>ncl 22>   printMinMax(u200,0)<br><b>(0)     min=-48.11269687686956   max=57.21568789939967</b><br>ncl 23>   print("------------------------------")<br>(0) ------------------------------<br>ncl 24>   printVarSummary(u850)<br><br>Variable: u850<br>Type: double<br>Total Size: 5504200 bytes<br>            688025 values<br>Number of Dimensions: 3<br>Dimensions and sizes: [time | 365] x [lat | 13] x [lon | 145]<br>Coordinates: <br>            time: [   0..524160]<br>            lat: [ -15..  15]<br>            lon: [   0.. 360]<br>Number Of Attributes: 1<br>  _FillValue :   -999000000<br>ncl 25>   printMinMax(u850,0)<br><b>(0)     min=-31.22578327383305   max=41.24650545117836</b></div><div><br></div><div>I want to run mjoclivar_14.ncl (<a href="https://www.ncl.ucar.edu/Applications/mjoclivar.shtml" target="_blank">https://www.ncl.ucar.edu/Applications/mjoclivar.shtml</a>), but I got is always both min and max data value -999000000. Could somebody help me making multivariate EOF with this file :</div><div><br></div><div>OLR anomaly data:<div class="gmail_chip gmail_drive_chip" style="width:396px;height:18px;max-height:18px;background-color:rgb(245,245,245);padding:5px;color:rgb(34,34,34);font-family:arial;font-style:normal;font-weight:bold;font-size:13px;border:1px solid rgb(221,221,221);line-height:1"><a href="https://drive.google.com/file/d/1bH67Flh4XgfYfSfsrBRPd0F7Gwq9qce5/view?usp=drive_web" style="display:inline-block;max-width:366px;overflow:hidden;text-overflow:ellipsis;white-space:nowrap;text-decoration:none;padding:1px 0px;border:medium none" target="_blank"><img style="vertical-align: bottom; border: medium none;" src="https://ssl.gstatic.com/docs/doclist/images/icon_10_generic_list.png"> <span dir="ltr" style="color:rgb(17,85,204);text-decoration:none;vertical-align:bottom">olranom.nc</span></a><img style="opacity: 0.55; float: right; display: none;"></div></div><div>200 hPa zonal wind anomaly data :<div class="gmail_chip gmail_drive_chip" style="width:396px;height:18px;max-height:18px;background-color:rgb(245,245,245);padding:5px;color:rgb(34,34,34);font-family:arial;font-style:normal;font-weight:bold;font-size:13px;border:1px solid rgb(221,221,221);line-height:1"><a href="https://drive.google.com/file/d/1Y_JmVgDatlxeQ2VsPSD4g-dh6Cc0XNna/view?usp=drive_web" style="display:inline-block;max-width:366px;overflow:hidden;text-overflow:ellipsis;white-space:nowrap;text-decoration:none;padding:1px 0px;border:medium none" target="_blank"><img style="vertical-align: bottom; border: medium none;" src="https://ssl.gstatic.com/docs/doclist/images/icon_10_generic_list.png"> <span dir="ltr" style="color:rgb(17,85,204);text-decoration:none;vertical-align:bottom">uanom.200.nc</span></a><img style="display: none; opacity: 0.55; float: right;"></div></div><div>850 hPa zonal wind anomaly data:<div class="gmail_chip gmail_drive_chip" style="width:396px;height:18px;max-height:18px;background-color:rgb(245,245,245);padding:5px;color:rgb(34,34,34);font-family:arial;font-style:normal;font-weight:bold;font-size:13px;border:1px solid rgb(221,221,221);line-height:1"><a href="https://drive.google.com/file/d/1TwaQL7zF8J3wIyCs8mmnQC1lqimPFsCV/view?usp=drive_web" style="display:inline-block;max-width:366px;overflow:hidden;text-overflow:ellipsis;white-space:nowrap;text-decoration:none;padding:1px 0px;border:medium none" target="_blank"><img style="vertical-align: bottom; border: medium none;" src="https://ssl.gstatic.com/docs/doclist/images/icon_10_generic_list.png"> <span dir="ltr" style="color:rgb(17,85,204);text-decoration:none;vertical-align:bottom">uanom.850.nc</span></a><img style="display: none; opacity: 0.55; float: right;"></div></div><div><br></div><div>This is my script for mjoclivar_14 :</div><div>ncl 0> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/diagnostics_cam.ncl"<br>dim_rmvmean_n( olr, 0)              ; (time,lon)<br>   u850  = dimncl 1> ;******************************************************<br>ncl 2> ;<br>ncl 3> ; mjoclivar_14.ncl<br>ncl 4> ;<br>ncl 5> ;***********************************************************<br>ncl 6> ; Combined EOFs<br>ncl 7> ; Latest Update: July, 2016: Eun-Pa Lim; Bureau of Meteorology, Australia<br>ncl 8> ;***********************************************************<br>ncl 9> ;;<br>ncl 10> ;;      The following are automatically loaded from 6.2.0 onward<br>ncl 11> ;;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"  <br>ncl 12> ;;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"  <br>ncl 13> ;;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" <br>ncl 14> <br>ncl 15> undef("read_rename")<br>ncl 16> function read_rename(f[1]:file, varName[1]:string       \<br>var_olr  = dim_avg_n_Wrap( var_olr , 0)     <br>  zavg_var_u850 = dim_avg_n_Wrap( var_u850, 0)<br>  zavg_var_u200 = dim_avg_n_Wrap(ncl 16>                     ,iStrt[1]:integer, iLast[1]:integer \<br>ncl 16>                     ,latS[1]:numeric , latN[1]:numeric  )<br>ncl 17> ; Utility to force specific named dimensions<br>ncl 18> ; This is done for historical reasons (convenience) <br>ncl 19> begin<br>ncl 20>    work    = f->$varName$(iStrt:iLast,{latS:latN},:)   ; (time,lat,lon)<br>ncl 21>    work!0  = "time"                                    ; CAM model names<br>ncl 22>    work!1  = "lat"<br>ncl 23>    work!2  = "lon"<br>ncl 24>    return(work)<br>ncl 25> end<br>ncl 26> ; =========================>  MAIN  <==============================<br>ncl 27> begin<br>ncl 28>    neof    =  2<br>ncl 29> <br>ncl 30>    latS    = -15<br>ncl 31>    latN    =  15<br>ncl 32> <br>mncl 33>    ymdStrt = 20180101                         ; start yyyymmdd<br>ncl 34>    ymdLast = 20181231                         ; last  <br>ncl 35> <br>yncl 36>    yrStrt  = ymdStrt/10000<br>ncl 37>    yrLast  = ymdLast/10000<br>ncl 38> <br>ancl 39>    pltDir  = "./"                             ; plot directory<br>ncl 40>    pltType = "png" <br>ncl 41>    pltName = "mjoclivar"   <br>ncl 42> <br>ncl 43>    diri    = "./"                             ; input directory<br>ncl 44> <br>oncl 45>    filolr  = "<a href="http://olranom.nc" target="_blank">olranom.nc</a>"<br>ncl 46>    filu200 = "<a href="http://uanom.200.nc" target="_blank">uanom.200.nc</a>"<br>ncl 47>    filu850 = "<a href="http://uanom.850.nc" target="_blank">uanom.850.nc</a>"<br>ncl 48> <br>ncl 49> ;************************************************<br>ncl 50> ; create BandPass Filter<br>ncl 51> ;************************************************<br>oncl 52>   ihp      = 2                             ; bpf=>band pass filter<br>ncl 53>   nWgt     = 201<br>ncl 54>   sigma    = 1.0                           ; Lanczos sigma<br>ncl 55>   fca      = 1./100.<br>ncl 56>   fcb      = 1./20.<br>ncl 57>   wgt      = filwgts_lanczos (nWgt, ihp, fca, fcb, sigma )<br>ncl 58> <br>*ncl 59> ;***********************************************************<br>ncl 60> ; Find the indices corresponding to the start/end times<br>ncl 61> ;***********************************************************<br>ncl 62>    f       = addfile (diri+filolr , "r")                         <br>ncl 63>    TIME    = f->time                          ; days since ...<br>ncl 64> <br>ncl 65>    YMD     = cd_calendar(TIME, -2)            ; entire (time,6)<br>ncl 66> <br>ncl 67>    iStrt   = ind(YMD.eq.ymdStrt)              ; index start<br>ncl 68>    iLast   = ind(YMD.eq.ymdLast)              ; index last <br>ncl 69>    delete([/ TIME, YMD /])<br>ncl 70> <br> ncl 71> ;***********************************************************<br>ncl 72> ; Read anomalies<br>ncl 73> ;***********************************************************<br>ncl 74> <br>encl 75>    work    = read_rename(f,"olranom",iStrt,iLast,latS,latN) ; (time,lat,lon) <br>ncl 76>    OLR     = dim_avg_n_Wrap(work, 1)                         ; (time,lon)<br>ncl 77>    delete(work)<br>ncl 78>    <br>ncl 79>    f       = addfile (diri+filu850 , "r")                         <br>ncl 80>    work    = read_rename(f,"uanom",iStrt,iLast,latS,latN) ; (time,lat,lon) <br>ncl 81>    U850    = dim_avg_n_Wrap(work, 1)          ; (time,lon)<br>ncl 82>    delete(work)<br>ncl 83>    <br>ncl 84>    f       = addfile (diri+filu200 , "r")                         <br>ncl 85>    work    = read_rename(f,"uanom",iStrt,iLast,latS,latN) ; (time,lat,lon) <br>ncl 86>    U200    = dim_avg_n_Wrap(work, 1)          ; (time,lon)<br>ncl 87> <br>0ncl 88>    dimw    = dimsizes( work )<br>ncl 89>    ntim    = dimw(0)<br>ncl 90>    nlat    = dimw(1)<br>)ncl 91>    mlon    = dimw(2)<br>uncl 92>    delete(work)<br>ncl 93> <br>oncl 94>    lon     = OLR&lon                                          <br>ncl 95>    time    = OLR&time                         <br>ncl 96>    date    = cd_calendar(time, -2)            ; yyyymmdd<br>ncl 97> <br>encl 98> ;************************************************<br>ncl 99> ; Apply the band pass filter to the original anomalies<br>ncl 100> ;************************************************<br>ncl 101>    olr   = wgt_runave_n_Wrap ( OLR, wgt, 0, 0) ; (time,lon)<br>ncl 102>    u850  = wgt_runave_n_Wrap (U850, wgt, 0, 0)<br>ncl 103>    u200  = wgt_runave_n_Wrap (U200, wgt, 0, 0)<br>ncl 104> <br>*ncl 105> ;************************************************<br>ncl 106> ; remove temporal means of band pass series: *not* necessary <br>ncl 107> ;************************************************<br>ncl 108>    olr   = dim_rmvmean_n( olr, 0)              ; (time,lon)<br>ncl 109>    u850  = dim_rmvmean_n(u850, 0)<br>ncl 110>    u200  = dim_rmvmean_n(u200, 0)<br>ncl 111> <br>nncl 112> ;************************************************<br>ncl 113> ; Compute the temporal variance at each lon<br>ncl 114> ;************************************************<br>ncl 115>    var_olr  = dim_variance_n_Wrap( olr, 0)     ; (lon)<br>ncl 116>    var_u850 = dim_variance_n_Wrap(u850, 0)<br>ncl 117>    var_u200 = dim_variance_n_Wrap(u200, 0)<br>ncl 118> <br>ncl 119> ;************************************************<br>ncl 120> ; Compute the zonal mean of the temporal variance<br>ncl 121> ;************************************************<br>ncl 122>   zavg_var_olr  = dim_avg_n_Wrap( var_olr , 0)     <br>ncl 123>   zavg_var_u850 = dim_avg_n_Wrap( var_u850, 0)<br>ncl 124>   zavg_var_u200 = dim_avg_n_Wrap( var_u200, 0)<br>ncl 125> <br>rncl 126> ;************************************************<br>ancl 127> ; Normalize by sqrt(avg_var*)<br>ncl 128> ;************************************************<br>ncl 129>   zsqrt = where(zavg_var_olr.gt.0,sqrt(zavg_var_olr),zavg_var_olr@_FillValue)<br>ncl 130>   olr   =  olr/zsqrt          ; (time,lon)<br>;ncl 131>   z850sqrt =where(zavg_var_u850.gt.0,sqrt(zavg_var_u850),zavg_var_u850@_FillValue)<br>ncl 132>   u850  = u850/z850sqrt<br>ncl 133>   z200sqrt= where(zavg_var_u200.gt.0,sqrt(zavg_var_u200),zavg_var_u200@_FillValue)<br>)ncl 134>   u200  = u200/z200sqrt<br>ncl 135> <br>cncl 136> ;************************************************<br>ncl 137> ; Combine the normalized data into one variable<br>ncl 138> ;************************************************<br>ncl 139>   cdata     = new ( (/3*mlon,ntim/), typeof(olr), getFillValue(olr))<br>ncl 140>   do ml=0,mlon-1<br>ncl 141>      cdata(ml       ,:) = (/  olr(:,ml) /)<br>ncl 142>      cdata(ml+  mlon,:) = (/ u850(:,ml) /)<br>ncl 143>      cdata(ml+2*mlon,:) = (/ u200(:,ml) /)<br>ncl 144>   end do<br>ncl 145> <br>=ncl 146> ;************************************************<br>ncl 147> ; Compute **combined** EOF; Sign of EOF is arbitrary<br>ncl 148> ;************************************************<br>ncl 149>   eof_cdata    = eofunc(cdata   , neof, False)      ; (neof,3*mlon)<br>ncl 150>   print("==============")<br>ncl 151>   printVarSummary(eof_cdata)<br>ncl 152>   printMinMax(eof_cdata, True)<br>ncl 153> <br>)ncl 154>   eof_ts_cdata = eofunc_ts(cdata,eof_cdata,False)   ; (neof,3*ntim)<br>ncl 155>   print("==============")                                  <br>ncl 156>   printVarSummary(eof_ts_cdata)<br>ncl 157>   printMinMax(eof_ts_cdata, True)<br>ncl 158> <br>ncl 159> ;************************************************<br>ncl 160> ; For clarity, explicitly extract each variable. Create time series <br>ncl 161> ;************************************************<br>ncl 162> <br>tncl 163>   nvar = 3  ; "olr", "u850", "u200"<br>ncl 164>   ceof = new( (/nvar,neof,mlon/), typeof(cdata), getFillValue(cdata))<br>ncl 165> <br>dncl 166>   do n=0,neof-1<br>ncl 167>      ceof(0,n,:) = eof_cdata(n,0:mlon-1)      ; olr<br>ncl 168>      ceof(1,n,:) = eof_cdata(n,mlon:2*mlon-1) ; u850<br>ncl 169>      ceof(2,n,:) = eof_cdata(n,2*mlon:)       ; u200<br>ncl 170>   end do<br>ncl 171> <br>(ncl 172>   ceof!0   = "var"<br>ncl 173>   ceof!1   = "eof"<br>ncl 174>   ceof!2   = "lon"   <br>ncl 175>   ceof&lon =  olr&lon<br>ncl 176> <br>incl 177>   ceof_ts        = new( (/nvar,neof,ntim/), typeof(cdata), getFillValue(cdata))<br>ncl 178>   ceof_ts(0,:,:) = eofunc_ts_Wrap( olr(lon|:,time|:),ceof(0,:,:),False)   ; (0,neof,ntim)<br>ncl 179>   ceof_ts(1,:,:) = eofunc_ts_Wrap(u850(lon|:,time|:),ceof(1,:,:),False)   ; (1,neof,ntim)<br>ncl 180>   ceof_ts(2,:,:) = eofunc_ts_Wrap(u200(lon|:,time|:),ceof(2,:,:),False)   ; (2,neof,ntim)<br>ncl 181> <br>tncl 182> ;**********************************************t*<br>ncl 183> ; Add code contributed by Marcus N. Morgan, Florida Institute of Technology; Feb 2015<br>ncl 184> ; Calculate % variance (pcv_ )accounted for by OLR, U850 and U200<br>ncl 185> ;************************************************<br>ncl 186> <br>ncl 187>     pcv_eof_olr  = new(neof,typeof(ceof))<br>ncl 188>     pcv_eof_u850 = new(neof,typeof(ceof))<br>ncl 189>     pcv_eof_u200 = new(neof,typeof(ceof))<br>ncl 190>       <br>ncl 191>     do n=0,neof-1<br>ncl 192>        pcv_eof_olr(n)  = avg((ceof(0,n,:)*sqrt(ceof@eval(n)))^2)*100<br>ncl 193>        pcv_eof_u850(n) = avg((ceof(1,n,:)*sqrt(ceof@eval(n)))^2)*100<br>ncl 194>        pcv_eof_u200(n) = avg((ceof(2,n,:)*sqrt(ceof@eval(n)))^2)*100<br>ncl 195>      ;;print("pcv: neof="+(n+1)+":  "+pcv_eof_olr(n)+"  "+pcv_eof_u850(n)+"  "+pcv_eof_u200(n))<br>ncl 196>     end do<br>ncl 197> <br>ncl 198> ;************************************************<br>ncl 199> ; Change sign of EOFs for spatial structures of PC1 and PC2 <br>ncl 200> ; to represent convection over the tropical Indian Ocean and the tropical western Pacific Ocean, respectively <br>ncl 201> ; (Ad hoc approach) <br>ncl 202> ;************************************************<br>ncl 203>                    <br>ncl 204>   imax_olr_eof1   = maxind(ceof(0,0,:))         <br>ncl 205>   imax_olr_eof2   = maxind(ceof(0,1,:))    <br>ncl 206> <br>)ncl 207>   lonmax_eof1 = ceof&lon(imax_olr_eof1)      ; longitude of max value (i.e. suppressed convection)<br>ncl 208>   lonmax_eof2 = ceof&lon(imax_olr_eof2)<br>ncl 209> <br>sncl 210>   if (lonmax_eof1.ge.60 .and. lonmax_eof1.lt.180) then  ; Change the sign of EOF1 <br>ncl 211>       ceof(:,0,:)       = -ceof(:,0,:)                  ; if OLR is positive<br>ncl 212>       ceof_ts(:,0,:)    = -ceof_ts(:,0,:)               ;  over the tropical Indian Ocean<br>ncl 213>       eof_cdata(0,:)    = -eof_cdata(0,:)<br></div><div>ncl 214>       eof_ts_cdata(0,:) = -eof_ts_cdata(0,:)<br>ncl 215>   end if<br>ncl 216> <br>encl 217>   if (lonmax_eof2.ge.120 .and. lonmax_eof2.lt.180) then  ; Change the sign of EOF2<br>ncl 218>       ceof(:,1,:)       = -ceof(:,1,:)                   ; if OLR is positive <br>ncl 219>       ceof_ts(:,1,:)    = -ceof_ts(:,1,:)                ; over the tropical western Pacific Ocean<br>ncl 220>       eof_cdata(1,:)    = -eof_cdata(1,:)<br>ncl 221>       eof_ts_cdata(1,:) = -eof_ts_cdata(1,:)<br>ncl 222>   end if<br>ncl 223> <br>ncl 224>   print("==============")<br>ncl 225>   printVarSummary(eof_cdata)<br>ncl 226>   printMinMax(eof_cdata, True)<br>nncl 227> <br>ncl 228> ;************************************************<br>ncl 229> ; Compute cross correlation of each variable's EOF time series at zero-lag<br>ncl 230> ;************************************************<br>ncl 231>   r_olr_u850  = escorc(ceof_ts(0,:,:) , ceof_ts(1,:,:) )  ; (neof)<br>ncl 232>   r_olr_u200  = escorc(ceof_ts(0,:,:) , ceof_ts(2,:,:) )<br>ncl 233>   r_u850_u200 = escorc(ceof_ts(1,:,:) , ceof_ts(2,:,:) )<br>ncl 234> <br>ncl 235>   print("==============")<br>ncl 236>   do n=0,neof-1<br>ncl 237>      print("neof="+n \<br>ncl 237>           +"  r_olr_u850=" +sprintf("%4.3f",r_olr_u850(n))  \<br>ncl 237>           +"  r_olr_u200=" +sprintf("%4.3f",r_olr_u200(n))  \<br>ncl 237>           +"  r_u850_u200="+sprintf("%4.3f",r_u850_u200(n)) )<br>ncl 238>   end do<br>ncl 239>   print("==============")<br>ncl 240> <br>ncl 241> ;************************************************<br>ncl 242> ; Compute cross correlation of the multivariate EOF; EOF 1 vs EOF 2<br>ncl 243> ;************************************************<br>ncl 244> <br>ncl 245>   mxlag     = 25<br>ncl 246>   rlag_01   = esccr(eof_ts_cdata(0,:),eof_ts_cdata(1,:), mxlag)   ; (N,mxlag+1)<br>ncl 247>   rlag_10   = esccr(eof_ts_cdata(1,:),eof_ts_cdata(0,:), mxlag)   ; (N,mxlag+1)<br>ncl 248>   ccr_12    = new ( (/2*mxlag+1/), float)    <br>ncl 249> <br>ncl 250>   ccr_12(mxlag:)    = rlag_10(0:mxlag)   <br>ncl 251>   ccr_12(0:mxlag)   = rlag_01(::-1)       ; reverse order<br>ncl 252> ;;print(ccr_12)<br>ncl 253> <br>ncl 254> <br>ncl 255> ;************************************************<br>ncl 256> ; Normalize the multivariate EOF 1&2 component time series<br>ncl 257> ; Compute (PC1^2+PC2^2): values > 1 indicate "strong" periods<br>encl 258> ;************************************************<br>ncl 259>   eof_ts_cdata(0,:) = eof_ts_cdata(0,:)/stddev(eof_ts_cdata(0,:))<br>ncl 260>   eof_ts_cdata(1,:) = eof_ts_cdata(1,:)/stddev(eof_ts_cdata(1,:))<br>ncl 261> <br>ncl 262>   mjo_ts_index      = eof_ts_cdata(0,:)^2 + eof_ts_cdata(1,:)^2 <br>ncl 263>   mjo_ts_index_smt  = runave(mjo_ts_index, 91, 0) ; 91-day running mean<br> ncl 264> <br>ncl 265>   nGood   = num(.not.ismissing(mjo_ts_index))     ; # non-missing<br>ncl 266>   nStrong = num(mjo_ts_index .ge. 1.0)<br>ncl 267>   print("nGood="+nGood+"   nStrong="+nStrong+"   nOther="+(nGood-nStrong))<br>ncl 268> <br>rncl 269> ;************************************************<br>uncl 270> ; Write PC results to netCDF for use in another example.<br>ncl 271> ;************************************************<br>ncl 272>   mjo_ts_index!0    = "time"<br>t<br>ncl 273>   mjo_ts_index&time = time <br>ncl 274>   mjo_ts_index@long_name = "MJO PC INDEX" <br><br><br>;****ncl 275>   mjo_ts_index@info      = "(PC1^2 + PC2^2)" <br>ncl 276> <br>ncl 277>   PC1           = eof_ts_cdata(0,:)<br>ncl 278>   PC1!0         = "time"<br>********************************************<br>; MJO "strong" index <br>;********************ncl 279>   PC1&time      =  time<br>ncl 280>   PC1@long_name = "PC1"<br>ncl 281>   PC1@info      = "PC1/stddev(PC1)"<br>ncl 282> <br>oncl 283>   PC2           = eof_ts_cdata(1,:)<br>ncl 284>   PC2!0         = "time"<br>ncl 285>   PC2&time      =  time<br>ncl 286>   PC2@long_name = "PC2"<br>ncl 287>   PC2@info      = "PC2/stddev(PC2)"<br>ncl 288> <br> ncl 289>   diro = "./"<br>ncl 290>   filo = "MJO_PC_INDEX.nc"<br>ncl 291>   system("/bin/rm -f "+diro+filo)   ; remove any pre-existing file<br>ncl 292>   ncdf = addfile(diro+filo,"c")     ; open output netCDF file<br>ncl 293>                                     ; make time an UNLIMITED dimension <br>ncl 294>   filedimdef(ncdf,"time",-1,True)   ; recommended  for most applications<br>ncl 295>                                     ; output variables directly<br>ncl 296>   ncdf->MJO_INDEX = mjo_ts_index    <br>ncl 297>   ncdf->PC1       = PC1     <br>ncl 298>   ncdf->PC2       = PC2     <br>ncl 299> <br>ncl 300> ;------------------------------------------------------------<br>ncl 301> ; PLOTS<br>ncl 302> ;------------------------------------------------------------<br>Pncl 303> <br>ncl 304>   yyyymmdd = cd_calendar(time, -2)<br>ncl 305>   yrfrac   = yyyymmdd_to_yyyyfrac(yyyymmdd, 0.0)<br>ncl 306>   delete([/ yrfrac@long_name, lon@long_name /])<br>ncl 307> <br>ncl 308>   day      = ispan(-mxlag, mxlag, 1)<br>ncl 309>  ;day@long_name = "lag (day)"<br>ncl 310> <br>ncl 311>   pltPath = pltDir+pltName<br>ncl 312> <br>ncl 313>   wks = gsn_open_wks(pltType,pltPath)<br>ncl 314>   gsn_define_colormap(wks,"default")<br>ncl 315>   plot = new(3,graphic)                <br>ncl 316> <br>ncl 317> ;************************************************<br>ncl 318> ; Multivariate EOF plots<br>ncl 319> ;************************************************<br>ncl 320>   rts           = True<br>ncl 321>   rts@gsnDraw   = False       ; don't draw yet<br>ncl 322>   rts@gsnFrame  = False       ; don't advance frame yet<br>ncl 323>   rts@gsnScale  = True        ; force text scaling               <br>ncl 324> <br>ncl 325>   rts@vpHeightF = 0.40        ; Changes the aspect ratio<br>ncl 326>   rts@vpWidthF  = 0.85<br>ncl 327>   rts@vpXF      = 0.10        ; change start locations<br>ncl 328>   rts@vpYF      = 0.75        ; the plot<br>ncl 329>   rts@xyLineThicknesses = (/2, 2, 2/)<br>ncl 330>   rts@xyLineColors      = (/"black","red","blue"/)<br>ncl 331>   rts@gsnYRefLine       = 0.                  ; reference line   <br>ncl 332>   rts@trXMaxF           = max(lon)<br>ncl 333>   rts@trXMinF           = min(lon)<br>ncl 334> <br>ncl 335>   rts@pmLegendDisplayMode    = "Always"            ; turn on legend<br>ncl 336>   rts@pmLegendSide           = "Top"               ; Change location of <br>ncl 337>   rts@pmLegendParallelPosF   = 1.16                ; move units right<br>ncl 338>   rts@pmLegendOrthogonalPosF = -0.50               ; move units down<br>ncl 339>   rts@pmLegendWidthF         = 0.15                ; Change width and<br>ncl 340>   rts@pmLegendHeightF        = 0.15                ; height of legend.<br>ncl 341>   rts@lgLabelFontHeightF     = 0.0175<br>ncl 342> <br>ncl 343> <br>ncl 344>   rtsP                       = True                ; modify the panel plot<br>ncl 345> ;  rtsP@gsnMaximize           = True                ; large format<br>ncl 346>   rtsP@gsnPanelMainString     = "Multivariate EOF: 15S-15N: "+yrStrt+"-"+yrLast <br>ncl 347>   <br>ncl 348>   do n=0,neof-1<br>ncl 349>     rts@xyExplicitLegendLabels = (/"OLR: "+sprintf("%4.1f", pcv_eof_u200(n)) +"%" \<br>ncl 349>                                   ,"U850: "+sprintf("%4.1f", pcv_eof_u850(n))+"%" \<br>ncl 349>                                   ,"U200: "+sprintf("%4.1f", pcv_eof_olr(n))+"%" /)<br>ncl 350>     rts@gsnLeftString  = "EOF "+(n+1)<br>ncl 351>     rts@gsnRightString = sprintf("%3.1f",ceof@pcvar(n))  +"%"<br>ncl 352>     plot(n) = gsn_csm_xy (wks,lon,ceof(:,n,:),rts)<br>ncl 353>   end do<br>ncl 354>   gsn_panel(wks,plot(0:1),(/2,1/),rtsP)     ; now draw as one plot<br>ncl 355> <br>ncl 356> ;-----------------------------------------<br>ncl 357> ; The following doesn't work with some older versions of NCL<br>ncl 358> ; With old versions, the user must delete each individually.<br>ncl 359> ;----------------------------------------- <br>ncl 360>   delete([/ rts@xyExplicitLegendLabels, rts@pmLegendDisplayMode \<br>ncl 360>           , rts@xyLineThicknesses     , rts@gsnLeftString       \<br>ncl 360>           , rts@gsnRightString        , rts@xyLineColors        \<br>ncl 360>           , rts@trXMaxF               , rts@trXMinF             /] )<br>ncl 361> <br>ncl 362>   lag                        = ispan(-mxlag,mxlag,1)<br>ncl 363>   lag@long_name              = "lag (days)"<br>ncl 364> <br>ncl 365>   plot(0)                    = gsn_csm_xy (wks, lag ,ccr_12,rts)<br>ncl 366>   rtsP@gsnPanelMainString    = "Cross Correlation: Multivariate EOF: 15S-15N: " \<br>ncl 366>                    +  yrStrt+"-"+yrLast <br>ncl 367>   rtsP@gsnPaperOrientation   = "portrait"        ; force portrait<br>ncl 368>   gsn_panel(wks,plot(0),(/1,1/),rtsP)     ; now draw as one plot<br>ncl 369> <br>ncl 370> ;************************************************<br>ncl 371> ; MJO "strong" index <br>ncl 372> ;************************************************<br>ncl 373>   rts@gsnYRefLine        = 1.0<br>ncl 374>   rts@gsnYRefLineColor   = "black"<br>ncl 375>   rts@xyMonoDashPattern  = True<br>ncl 376>   rts@xyLineColors       = (/"black", "blue"/)<br>ncl 377>   rts@xyLineThicknesses  = (/1, 2/)<br>ncl 378>   rts@pmLegendDisplayMode    = "Always"            ; turn on legend<br>ncl 379>   rts@pmLegendWidthF         = 0.12                ; Change width and<br>ncl 380>   rts@pmLegendHeightF        = 0.10                ; height of legend.<br>ncl 381>   rts@pmLegendParallelPosF   = 0.86                ; move units right<br>ncl 382>   rts@pmLegendOrthogonalPosF = -0.40               ; move units down<br>ncl 383>   rts@xyExplicitLegendLabels = (/"daily", "91-day runavg" /)<br>ncl 384> <br>ncl 385>   mjo_ind_plt = new ( (/2,ntim/), typeof(mjo_ts_index))<br>ncl 386>   mjo_ind_plt(0,:) = mjo_ts_index<br>ncl 387>   mjo_ind_plt(1,:) = (/ mjo_ts_index_smt /)<br>ncl 388>   plot(0) = gsn_csm_xy(wks, yrfrac,mjo_ind_plt,rts)<br>ncl 389> <br>ncl 390>   rtsP@gsnPanelMainString   = "MJO Index: (PC1^2+ PC2^2) : 15S-15N: "+yrStrt+"-"+yrLast <br>ncl 391>   gsn_panel(wks,plot(0),(/1,1/),rtsP)     ; now draw as one plot<br>ncl 392> <br>ncl 393>  end<br><br><br> sub drveof: ier,jopt=  10  0 returned from vcmssm/crmssm<br>warning:eofunc: 10 eigenvectors failed to converge<br>(0)     ==============<br><br>Variable: eof_cdata<br>Type: double<br>Total Size: 6960 bytes<br>            870 values<br>Number of Dimensions: 2<br>Dimensions and sizes:       [2] x [435]<br>Coordinates: <br>Number Of Attributes: 5<br>  _FillValue :    -999000000<br>  method : no transpose<br>  matrix :       covariance<br>  pcvar :  ( -9.99e+08, -9.99e+08 )<br>  eval :     ( -999000000, -999000000 )<br>(0)  <br>(0)  min=-999000000   max=-999000000<br>warning:eofunc_ts: 2 eigenvectors failed to converge<br>(0)     ==============<br><br>Variable: eof_ts_cdata<br>Type: double<br>Total Size: 5840 bytes<br>            730 values<br>Number of Dimensions: 2<br>Dimensions and sizes:    [2] x [365]<br>Coordinates: <br>Number Of Attributes: 3<br>  _FillValue :    -999000000<br>  matrix : covariance<br>  ts_mean :        ( -999000000, -999000000 )<br>(0)  <br>(0)  min=-999000000   max=-999000000<br><b>warning:eofunc_ts: 2 eigenvectors failed to converge<br>warning:eofunc_ts: 2 eigenvectors failed to converge<br>warning:eofunc_ts: 2 eigenvectors failed to converge<br>warning:avg: Entire input array contains missing values, can't compute average<br>warning:avg: Entire input array contains missing values, can't compute average<br>warning:avg: Entire input array contains missing values, can't compute average<br>warning:avg: Entire input array contains missing values, can't compute average<br>warning:avg: Entire input array contains missing values, can't compute average<br>warning:avg: Entire input array contains missing values, can't compute average<br>fatal:Subscript out of range, error in subscript #0<br>fatal:An error occurred reading unknown<br>fatal:["Execute.c":8635]:Execute: Error occurred at or near line 207</b></div><div><br></div><div>Thank you. <br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Fri, Dec 20, 2019 at 1:15 PM Barry Lynn <<a href="mailto:barry.h.lynn@gmail.com" target="_blank">barry.h.lynn@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr">Hi:<div><br></div><div>Could you please define what you mean by "missing."  Did you PrntMinMax(XXXX,False)?</div><div><br></div><div><br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Fri, Dec 20, 2019 at 1:20 AM Rahpeni Fajarianti via ncl-talk <<a href="mailto:ncl-talk@ucar.edu" target="_blank">ncl-talk@ucar.edu</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>Hi Dennis, Thank you for your time. <br></div><div><br></div><div>I am going to get OLR and (850 hPa and 200 hPa) zonal wind anomaly  for year 2018 using mjoclivar_2.ncl (<a href="https://www.ncl.ucar.edu/Applications/Scripts/mjoclivar_2.ncl" target="_blank">https://www.ncl.ucar.edu/Applications/Scripts/mjoclivar_2.ncl</a>) to create netcdf file as input data for mjoclivar_14. <br></div><div><br></div><div>The dataset I used were OLR daily mean (<a href="https://www.esrl.noaa.gov/psd/data/gridded/data.interp_OLR.html" target="_blank">https://www.esrl.noaa.gov/psd/data/gridded/data.interp_OLR.html</a>) and u wind daily mean (<a href="https://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.pressure.html" target="_blank">https://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.pressure.html</a>). <br></div><div><br></div><div>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.</div><div>can anybody help to fix this? Thank you<br></div><div><br></div><div><b>This is what I got for OLR anomaly:</b></div><div>Variable: f2<br>Type: file<br>filename: olr.day.anomalies.2018<br>path:   /home/peni/<a href="http://olr.day.anomalies.2018.nc" target="_blank">olr.day.anomalies.2018.nc</a><br>   file global attributes:<br>      creation_date : Sel Des 17 08:24:41 WIB 2019<br>      source_file : <a href="http://olr.day.mean.nc" target="_blank">olr.day.mean.nc</a><br>      title : OLR: Daily Anomalies: 2018-2018<br>   dimensions:<br>      time = 365  // unlimited<br>      lat = 73<br>      lon = 144<br>   variables:<br>      double time ( time )<br>         units :     hours since 1800-01-01 00:00:0.0<br>         long_name :     Time<br>         actual_range :      ( 1528872, 1919688 )<br>         delta_t :   0000-00-01 00:00:00<br>         standard_name :      time<br>         axis :      T<br><br>      float lat ( lat )<br>         units :  degrees_north<br>         actual_range :     ( 90, -90 )<br>         long_name :  Latitude<br>         standard_name : latitude<br>         axis :  Y<br><br>      float lon ( lon )<br>         units :  degrees_east<br>         long_name : Longitude<br>         actual_range : (  0, 360 )<br>         standard_name :     longitude<br>         axis : X<br><br>      float OLR_anom ( time, lat, lon )<br>         _FillValue :     32766<br>         long_name :        Anomalies: Daily OLR<br>         units :     W/m^2<br><br>      float OLR_anom_sm ( time, lat, lon )<br>         _FillValue :      32766<br>         long_name :        Anomalies from Smooth Daily Climatology<br>         units :  W/m^2</div><div><br></div><div><br></div><div><b>this is what I got for zonal wind anomaly:</b></div><div>Variable: f3<br>Type: file<br>filename:   uwnd.anomalies<br>path:   /home/peni/<a href="http://uwnd.anomalies.nc" target="_blank">uwnd.anomalies.nc</a><br>   file global attributes:<br>      creation_date : Sel Des 17 12:56:50 WIB 2019<br>      source_file : <a href="http://uwnd.day.mean.200.nc" target="_blank">uwnd.day.mean.200.nc</a><br>   dimensions:<br>      time = 365  // unlimited<br>      lat = 73<br>      lon = 145<br>   variables:<br>      double time ( time )<br>         long_name : Time<br>         units :     minutes since 2018-01-01 00:00<br><br>      double lat ( lat )<br>         units :    degrees_north<br>         long_name :        Latitude<br><br>      double lon ( lon )<br>         units :  degrees_east<br>         long_name : Longitude<br><br>      float U_anom ( time, lat, lon )<br>         long_name :        Anomalies from Daily Climatology<br>         _FillValue :    -9.99e+08<br><br>      float U_anom_sm ( time, lat, lon )<br>         long_name :     Anomalies from Smooth Daily Climatology<br>         _FillValue :     -9.99e+08</div><div><br></div><div><br></div><div><b>This is my script:</b></div><div>;-------------------------------------------------------------<br>; mjoclivar_2.ncl<br>;-------------------------------------------------------------<br>; These files are loaded by default in NCL V6.2.0 and newer<br>; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"  <br>; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"  <br>; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" <br>;******************************************************<br><br>begin<br>;******************************************************<br>; User specifications<br>;******************************************************<br>   NC      = True                             ; create netCDF?       <br>   PLOT    = False                            ; sample plots?<br><br>   ymdStrt = 20180101                         ; start yyyymmdd<br>   ymdLast = 20181231                         ; last  <br>   yrStrt  = ymdStrt/10000<br>   yrLast  = ymdLast/10000<br>  <br>   nhar    = 4                                ; number of fourier comp      <br>   <br>   var     = "uwnd2"                            ; name of file                   <br>   vName   = "U"                               ; name for plots<br><br>  ;diri    = "/project/cas/shea/CDC/"         ; input directory   <br>  ;diri    = "/Users/shea/Data/AMWG/"         ; input directory   <br>   diri    = "/home/peni/"                    ; new input directory<br>   fili    = var+".<a href="http://day.mean.850.nc" target="_blank">day.mean.850.nc</a>"               ; input file       <br><br>   if (NC) then<br>       diro= "/home/peni/"         ; output dir<br>      ;filo= var+".<a href="http://day.anomalies.nc" target="_blank">day.anomalies.nc</a>"          ; output file  <br>       filo= var+".<a href="http://anomalies.nc" target="_blank">anomalies.nc</a>"  ; output file  <br>   end if<br><br>   if (PLOT) then<br>       wksType = "png"                        ; send graphics to PNG file<br>       wksName = "mjoclivar"                  ; "mjo."+yrStrt+"_"+yrLast<br>   end if<br><br>;***********************************************************<br>; Read user specified time and create required yyyyddd                    <br>;***********************************************************<br>   f       = addfile (diri+fili , "r")                          <br><br>   time    = f->time                          ; all times on file<br>   ymd     = cd_calendar(time, -2)            ; yyyymmdd<br>   iStrt   = ind(ymd.eq.ymdStrt)              ; index start<br>   iLast   = ind(ymd.eq.ymdLast)              ; index last <br>   delete(time)<br>   delete(ymd)<br><br>;***********************************************************<br>; Read user specified time and create required yyyyddd                    <br>;***********************************************************<br>   time    = f->time(iStrt:iLast)             ; time:units = "hours since"<br>   TIME    = cd_calendar(time, 0)             ; type float <br>   year    = floattointeger( TIME(:,0) )<br>   month   = floattointeger( TIME(:,1) )<br>   day     = floattointeger( TIME(:,2) ) <br>   ddd     = day_of_year(year, month, day) <br>   yyyyddd = year*1000 + ddd                  ; needed for input<br><br>;***********************************************************<br>; Read data: dble2flt                                     <br>;*********************************************************** <br>   x       =  dble2flt( f->$var$(iStrt:iLast,:,:) )    ; convert to float <br>   printVarSummary( x ) <br><br>;***********************************************************<br>; Compute daily climatology: raw and then 'smoothed'  <br>;***********************************************************<br>   xClmDay = clmDayTLL(x, yyyyddd)     ; daily climatology at each grid point<br>   printVarSummary(xClmDay)            <br><br>;***********************************************************<br>; Compute smoothed daily climatology using 'nhar' harmonics<br>;***********************************************************<br>   xClmDay_sm = smthClmDayTLL(xClmDay, nhar)<br>   printVarSummary(xClmDay_sm)<br><br>;***********************************************************<br>; Compute daily anomalies using raw and smoothed climatologies<br>;***********************************************************<br>    xAnom      = calcDayAnomTLL (x, yyyyddd, xClmDay)     <br>    printVarSummary(xAnom)<br>    printMinMax(xAnom, True)<br><br>    xAnom_sm   = calcDayAnomTLL (x, yyyyddd, xClmDay_sm)     <br>    xAnom_sm@long_name = "Anomalies from Smooth Daily Climatology"<br>    printVarSummary(xAnom_sm)<br>    printMinMax(xAnom_sm, True)<br><br>    delete( x )    ; no longer needed<br><br><br>;***********************************************************<br>; Create netCDF: convenience use 'simple' method<br>;***********************************************************<br><br>    dimx   = dimsizes(xAnom)<br>    ntim   = dimx(0)<br>    nlat   = dimx(1)<br>    mlon   = dimx(2)<br><br>    if (NC) then<br><br>        lat    = f->lat<br>        lon    = f->lon<br><br>        system("/bin/rm -f "+diro+filo)      ; rm any pre-exist file, if any<br>        fnc    = addfile (diro+filo, "c")<br><br>        filAtt = 0<br>        filAtt@title         = vName+": Daily Anomalies: "+yrStrt+"  <br>     filAtt@source_file   = fili<br>        filAtt@creation_date = systemfunc("date")<br>        fileattdef( fnc, filAtt )         ; copy file attributes  <br><br>        setfileoption(fnc,"DefineMode",True)<br>        <br>        varNC      = vName+"_anom"<br>        varNC_sm   = vName+"_anom_sm"<br><br>        dimNames = (/"time", "lat", "lon"/)  <br> dimSizes = (/ -1   ,  nlat,  mlon/) <br>       dimUnlim = (/ True , False, False/)   <br>       filedimdef(fnc,dimNames,dimSizes,dimUnlim)<br><br>        filevardef(fnc, "time"  ,typeof(time),getvardims(time)) <br>        filevardef(fnc, "lat"   ,typeof(lat) ,getvardims(lat)) <br>        filevardef(fnc, "lon"   ,typeof(lon) ,getvardims(lon))<br>        filevardef(fnc, varNC   ,typeof(xAnom)   ,getvardims(xAnom))    <br>        filevardef(fnc, varNC_sm,typeof(xAnom)   ,getvardims(xAnom))    <br><br>        filevarattdef(fnc,"time"  ,time)                    ; copy time attributes<br>        filevarattdef(fnc,"lat"   ,lat)                     ; copy lat attributes<br>        filevarattdef(fnc,"lon"   ,lon)                     ; copy lon attributes<br>        filevarattdef(fnc,varNC   ,xAnom)                <br>        filevarattdef(fnc,varNC_sm,xAnom_sm)                <br><br>        fnc->time       = (/time/)     <br>        fnc->lat        = (/lat/)<br>        fnc->lon        = (/lon/)<br>        fnc->$varNC$    = (/xAnom   /)<br>        fnc->$varNC_sm$ = (/xAnom_sm/)<br>    end if</div><div>end</div><div><br></div><div><br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Thu, Dec 19, 2019 at 9:59 AM Rahpeni Fajarianti <<a href="mailto:rahpenifajarianti@gmail.com" target="_blank">rahpenifajarianti@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="auto"></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">---------- Forwarded message ---------<br>Dari: <strong class="gmail_sendername" dir="auto">Rahpeni Fajarianti</strong> <span dir="auto"><<a href="mailto:rahpenifajarianti@gmail.com" target="_blank">rahpenifajarianti@gmail.com</a>></span><br>Date: Kam, 19 Des 2019 9:10 AM<br>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<br>To: reza tisa <<a href="mailto:rezatisa1@gmail.com" target="_blank">rezatisa1@gmail.com</a>><br></div><br><br><div dir="auto"></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">---------- Forwarded message ---------<br>Dari: <strong class="gmail_sendername" dir="auto">Dennis Shea</strong> <span dir="auto"><<a href="mailto:shea@ucar.edu" rel="noreferrer" target="_blank">shea@ucar.edu</a>></span><br>Date: Kam, 19 Des 2019 7:57 AM<br>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<br>To: Rahpeni Fajarianti <<a href="mailto:rahpenifajarianti@gmail.com" rel="noreferrer" target="_blank">rahpenifajarianti@gmail.com</a>><br>Cc: Ncl-talk <<a href="mailto:ncl-talk@ucar.edu" rel="noreferrer" target="_blank">ncl-talk@ucar.edu</a>><br></div><br><br><div dir="ltr"><div>The GOLDEN RULE of data processing is: <b><span style="color:rgb(255,0,0)">***LOOK AT YOUR DATA***</span></b></div><div><br></div><div>This is *user* responsibility.  <b>*ALL* the data are 0.0</b></div><div><br></div><div><b>%></b> ncl rdData.ncl</div><div>=================</div><div>Variable: olr<br>Dimensions and sizes:  [time | 365] x [lat | 13] x [lon | 144]<br>Coordinates: <br>            time: [1910952..1919688]<br>            lat: [-15..15]<br>            lon: [ 0..357.5]<br>(0)   A<span style="color:rgb(0,0,255)"><b>nomalies: Daily OLR (W/m^2) : min=0   max=0</b></span><br>(0)     ------------------------------<br><br>Variable: u200<br>Dimensions and sizes: [time | 365] x [lat | 13] x [lon | 145]<br>Coordinates: <br>            time: [   0..524160]<br>            lat: [ -15..  15]<br>            lon: [   0.. 360]<br>Number Of Attributes: 2<br>  long_name :    Anomalies from Daily Climatology<br>  _FillValue :       -9.99e+08<br>(0)  <span style="color:rgb(0,0,255)"><b>Anomalies from Daily Climatology : min=0   max=0</b></span><br>(0) ------------------------------<br><br>Variable: u850<br>Dimensions and sizes: [time | 365] x [lat | 13] x [lon | 145]<br>Coordinates: <br>            time: [   0..524160]<br>            lat: [ -15..  15]<br>            lon: [   0.. 360]<br>Number Of Attributes: 2<br>  _FillValue :   -9.99e+08<br> <br>(0)      <b><span style="color:rgb(0,0,255)">Anomalies from Daily Climatology : min=0   max=0</span></b><br></div><div><br></div><br><div><br></div><div><br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Wed, Dec 18, 2019 at 4:58 PM Rahpeni Fajarianti <<a href="mailto:rahpenifajarianti@gmail.com" rel="noreferrer noreferrer" target="_blank">rahpenifajarianti@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>Thank you so much for responding. I attach some files here. The data and script that I used. <br></div><div>The data : <a href="http://olr.day.anomalies.2018.nc" rel="noreferrer noreferrer" target="_blank">olr.day.anomalies.2018.nc</a> , <a href="http://uwnd.anomalies.nc" rel="noreferrer noreferrer" target="_blank">uwnd.anomalies.nc</a> (for 200hPa u wind anomalies) , <a href="http://uwnd2.anomalies.nc" rel="noreferrer noreferrer" target="_blank">uwnd2.anomalies.nc</a> (for 850hPa u wind anomalies)</div><div>my script: mjoclivar_14.ncl</div><div><br></div><div>I am very new in NCL and very pleased for help. Thank you.<br></div><div><div class="gmail_chip gmail_drive_chip" style="width:396px;height:18px;max-height:18px;background-color:rgb(245,245,245);padding:5px;color:rgb(34,34,34);font-family:arial;font-style:normal;font-weight:bold;font-size:13px;border:1px solid rgb(221,221,221);line-height:1"><a href="https://drive.google.com/file/d/1nstwpHwEVEB9nbuWQPnB9nntN2YVBz79/view?usp=drive_web" style="display:inline-block;max-width:366px;overflow:hidden;text-overflow:ellipsis;white-space:nowrap;text-decoration:none;padding:1px 0px;border:medium none" rel="noreferrer noreferrer" target="_blank"><img style="vertical-align: bottom; border: medium none;" src="https://ssl.gstatic.com/docs/doclist/images/icon_10_generic_list.png"> <span dir="ltr" style="color:rgb(17,85,204);text-decoration:none;vertical-align:bottom">mjoclivar_14.ncl</span></a><img style="opacity: 0.55; float: right; display: none;"></div><div class="gmail_chip gmail_drive_chip" style="width:396px;height:18px;max-height:18px;background-color:rgb(245,245,245);padding:5px;color:rgb(34,34,34);font-family:arial;font-style:normal;font-weight:bold;font-size:13px;border:1px solid rgb(221,221,221);line-height:1"><a href="https://drive.google.com/file/d/1RQIk85UqimEzDzqWZ9V97ZpbqRG8h-cz/view?usp=drive_web" style="display:inline-block;max-width:366px;overflow:hidden;text-overflow:ellipsis;white-space:nowrap;text-decoration:none;padding:1px 0px;border:medium none" rel="noreferrer noreferrer" target="_blank"><img style="vertical-align: bottom; border: medium none;" src="https://ssl.gstatic.com/docs/doclist/images/icon_10_generic_list.png"> <span dir="ltr" style="color:rgb(17,85,204);text-decoration:none;vertical-align:bottom">olr.day.anomalies.2018.nc</span></a><img style="display: none; opacity: 0.55; float: right;"></div><div class="gmail_chip gmail_drive_chip" style="width:396px;height:18px;max-height:18px;background-color:rgb(245,245,245);padding:5px;color:rgb(34,34,34);font-family:arial;font-style:normal;font-weight:bold;font-size:13px;border:1px solid rgb(221,221,221);line-height:1"><a href="https://drive.google.com/file/d/1ZxF2tZnCr7q1LFQWtGhMOLvaIE3mwK-p/view?usp=drive_web" style="display:inline-block;max-width:366px;overflow:hidden;text-overflow:ellipsis;white-space:nowrap;text-decoration:none;padding:1px 0px;border:medium none" rel="noreferrer noreferrer" target="_blank"><img style="vertical-align: bottom; border: medium none;" src="https://ssl.gstatic.com/docs/doclist/images/icon_10_generic_list.png"> <span dir="ltr" style="color:rgb(17,85,204);text-decoration:none;vertical-align:bottom">uwnd2.anomalies.nc</span></a><img style="display: none; opacity: 0.55; float: right;"></div><div class="gmail_chip gmail_drive_chip" style="width:396px;height:18px;max-height:18px;background-color:rgb(245,245,245);padding:5px;color:rgb(34,34,34);font-family:arial;font-style:normal;font-weight:bold;font-size:13px;border:1px solid rgb(221,221,221);line-height:1"><a href="https://drive.google.com/file/d/1rm-f84yy1PkKlA8ZBG87nRGDbnD8auqx/view?usp=drive_web" style="display:inline-block;max-width:366px;overflow:hidden;text-overflow:ellipsis;white-space:nowrap;text-decoration:none;padding:1px 0px;border:medium none" rel="noreferrer noreferrer" target="_blank"><img style="vertical-align: bottom; border: medium none;" src="https://ssl.gstatic.com/docs/doclist/images/icon_10_generic_list.png"> <span dir="ltr" style="color:rgb(17,85,204);text-decoration:none;vertical-align:bottom">uwnd.anomalies.nc</span></a><img style="display: none; opacity: 0.55; float: right;"></div></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Thu, Dec 19, 2019 at 2:22 AM Dennis Shea <<a href="mailto:shea@ucar.edu" rel="noreferrer noreferrer" target="_blank">shea@ucar.edu</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>This is offlineI</div><div><br></div><div><div>If you do not know NCL well, please read the excellent tutorial at <br></div><div><a href="http://www.ncl.ucar.edu/Document/Manuals/NCL_User_Guide/" rel="noreferrer noreferrer" target="_blank"><b>http://www.ncl.ucar.edu/Document/Manuals/NCL_User_Guide/</b></a></div><div>---</div></div><div> 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.<br></div><div>e<br></div><div>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.</div><div>---<br></div><div><div>This is the result:</div><div>Variable: olr<br>Type: float<br></div><div>[SNIP]<br></div><div>Number Of Attributes: 5<br>  _FillValue :    32766<br><br></div><div>(0,0)   32766    ,==== All values are missing [_FillValue ]<br></div><div>[SNIP]    <br></div><div><br></div><div>===</div><div>Well .... at what point in your script did things go 'wrong'?</div><div>Use printVarSummary(), printMinMax(), etc</div><div>===</div><div><br></div><div>I will look at the issue(s) if you make your data available via (say) google drive or ftp.</div><div><br></div><div>ftp <a href="http://ftp.cgd.ucar.edu" rel="noreferrer noreferrer" target="_blank">ftp.cgd.ucar.edu</a></div><div>anonymous</div><div>your_email</div><div>cd incoming</div><div>put ...</div><div>put ...</div><div>put ...</div><div>quit<br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Wed, Dec 18, 2019 at 8:14 AM Rahpeni Fajarianti via ncl-talk <<a href="mailto:ncl-talk@ucar.edu" rel="noreferrer noreferrer" target="_blank">ncl-talk@ucar.edu</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>Thanks for your help. I had use print and printVarSummary statement. <br></div><div><br></div><div>This is the result:</div><div>Variable: olr<br>Type: float<br>Total Size: 210240 bytes<br>            52560 values<br>Number of Dimensions: 2<br>Dimensions and sizes:       [time | 365] x [lon | 144]<br>Coordinates: <br>            time: [1910952..1919688]<br>            lon: [ 0..357.5]<br>Number Of Attributes: 5<br>  _FillValue :     32766<br>  units :       W/m^2<br>  long_name :   Anomalies: Daily OLR<br>  average_op_ncl :       dim_avg_n over dimension(s): lat<br>  wgt_runave_op_ncl :        wgt_runave_n<br>(0,0)     32766<br>(0,1)    32766<br>(0,2)    32766<br>(0,3)    32766<br>(0,4)    32766<br>(0,5)    32766<br>(0,6)    32766<br>(0,7)    32766<br>(0,8)    32766<br>(0,9)    32766<br>(0,10)   32766<br>(0,11)   32766<br>(0,12)   32766<br>(0,13)   32766<br>(0,14)   32766<br>(0,15)   32766<br>(0,16)   32766<br>(0,17)   32766<br>(0,18)   32766<br>(0,19)   32766<br>(0,20)   32766<br>(0,21)   32766<br>(0,22)   32766<br>(0,23)   32766<br>(0,24)   32766<br>(0,25)   32766<br>(0,26)   32766<br>(0,27)   32766<br>(0,28)   32766<br>(0,29)   32766<br>(0,30)   32766<br>(0,31)   32766<br>(0,32)   32766<br>(0,33)   32766<br>(0,34)   32766<br>(0,35)   32766<br>(0,36)   32766<br>(0,37)   32766<br>(0,38)   32766<br>(0,39)   32766<br>(0,40)   32766<br>(0,41)   32766<br>(0,42)   32766<br>(0,43)   32766<br>(0,44)   32766<br>(0,45)   32766<br>(0,46)   32766<br>(0,47)   32766<br>(0,48)   32766<br>(0,49)   32766<br>(0,50)   32766<br>(0,51)   32766<br>(0,52)   32766<br>(0,53)   32766<br>(0,54)   32766<br>(0,55)   32766<br>(0,56)   32766<br>(0,57)   32766<br>(0,58)   32766<br>(0,59)   32766<br>(0,60)   32766<br>(0,61)   32766<br>(0,62)   32766<br>(0,63)   32766<br>(0,64)   32766<br>(0,65)   32766<br>(0,66)   32766<br>(0,67)   32766<br>(0,68)   32766<br>(0,69)   32766<br>(0,70)   32766<br>(0,71)   32766<br>(0,72)   32766<br>(0,73)   32766<br>(0,74)   32766<br>(0,75)   32766<br>(0,76)   32766<br>(0,77)   32766<br>(0,78)   32766<br>(0,79)   32766<br>(0,80)   32766<br>(0,81)   32766<br>(0,82)   32766<br>(0,83)   32766<br>(0,84)   32766<br>(0,85)   32766<br>(0,86)   32766<br>(0,87)   32766<br>(0,88)   32766<br>(0,89)   32766<br>(0,90)   32766<br>(0,91)   32766<br>(0,92)   32766<br>(0,93)   32766<br>(0,94)   32766<br>(0,95)   32766<br>(0,96)   32766<br>(0,97)   32766<br>(0,98)   32766<br>(0,99)   32766<br>(0,100)  32766<br>(0,101)  32766<br>(0,102)  32766<br>(0,103)  32766<br>(0,104)  32766<br>(0,105)  32766<br>(0,106)  32766<br>(0,107)  32766<br>(0,108)  32766<br>(0,109)  32766<br>(0,110)  32766<br>(0,111)  32766<br>(0,112)  32766<br>(0,113)  32766<br>(0,114)  32766<br>(0,115)  32766<br>(0,116)  32766<br>(0,117)  32766<br>(0,118)  32766<br>(0,119)  32766<br>(0,120)  32766<br>(0,121)  32766<br>(0,122)  32766<br>(0,123)  32766<br>(0,124)  32766<br>(0,125)  32766<br>(0,126)  32766<br>(0,127)  32766<br>(0,128)  32766<br>(0,129)  32766<br>(0,130)  32766<br>(0,131)  32766<br>(0,132)  32766<br>(0,133)  32766<br>(0,134)  32766<br>(0,135)  32766<br>(0,136)  32766<br>(0,137)  32766<br></div><div>(0,138)   32766 ect..</div><div><br></div><div>What should I fix to my script? <br></div><div>Thanks. <br></div><div><br></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Wed, Dec 18, 2019 at 8:50 AM Adam Phillips <<a href="mailto:asphilli@ucar.edu" rel="noreferrer noreferrer" target="_blank">asphilli@ucar.edu</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr">Hi Rahpeni,<div>Please always include ncl-talk on your replies; that way others can assist and see the correspondence. </div><div><br></div><div>You stated:</div><div><div>Thank you so much, it helps ! </div><div><br></div><div>I am still getting some errors.<br></div><div>I got :</div><div>fatal:Subscript out of range, error in subscript #1<br>fatal:An error occurred reading olr<br>fatal:["Execute.c":8635]:Execute: Error occurred at or near line 141</div></div><div><br></div><div><br></div><div>Note that it is for you to attempt to diagnose the issue before responding to ncl-talk. The error message is telling you what the problem is: You are attempting to access an index of a dimension that is bigger than the dimension itself. The lines around 141 are (from the email you sent me):</div><div><br></div><div>ncl 140>   do ml=0,mlon-1<br>ncl 141>      cdata(ml       ,:) = (/  olr(:,ml) /)<br>ncl 142>      cdata(ml+  mlon,:) = (/ u850(:,ml) /)<br>ncl 143>      cdata(ml+2*mlon,:) = (/ u200(:,ml) /)<br>ncl 144>   end do<br></div><div><br></div><div>NCL is zero based, so error in subscript #1 is referring to the 2nd index, line 141 shows that you are accessing the 2nd index of olr like this: (/  olr(:,ml) /)</div><div><br></div><div>ml at some point is bigger than the 2nd dimension size of olr. You will have to figure out why. Use print and printVarSummary statements to help to diagnose the issue. </div><div><br></div><div>Hope that helps!</div><div>Adam</div><div><br></div><div><br></div><div><br></div><div><br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Tue, Dec 17, 2019 at 1:40 AM Rahpeni Fajarianti via ncl-talk <<a href="mailto:ncl-talk@ucar.edu" rel="noreferrer noreferrer" target="_blank">ncl-talk@ucar.edu</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>Hi all ! <br></div><div>I want to make MJO clivar 14, but I have some problem.<br></div><div>This is :</div><div><div>fatal:divide: Division by 0, Can't continue</div>fatal:Div: operator failed, can't continue<br>fatal:["Execute.c":8635]:Execute: Error occurred at or near line 129</div><div><br></div><div>This is my script:<br></div><div>ncl 0> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/diagnostics_cam.ncl"</div>ncl 1> ;******************************************************<br> u850  = dim_rmvmean_n(u850, 0)<br>   u200  = dim_rmvmean_nncl 2> ;<br>(uncl 3> ; mjoclivar_14.ncl<br>ncl 4> ;<br>ncl 5> ;***********************************************************<br>ncl 6> ; Combined EOFs<br>hncl 7> ; Latest Update: July, 2016: Eun-Pa Lim; Bureau of Meteorology, Australia<br>ncl 8> ;***********************************************************<br>iancncl 9> ;;<br>nncl 10> ;;      The following are automatically loaded from 6.2.0 onward<br>ncl 11> ;;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"  <br>ncl 12> ;;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"  <br>ncl 13> ;;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" <br>ncl 14> <br>ncl 15> undef("read_rename")<br>ncl 16> function read_rename(f[1]:file, varName[1]:string       \<br>ncl 16>                     ,iStrt[1]:integer, iLast[1]:integer \<br>ncl 16>                     ,latS[1]:numeric , latN[1]:numeric  )<br>ncl 17> ; Utility to force specific named dimensions<br>ncl 18> ; This is done for historical reasons (convenience) <br>ncl 19> begin<br>ncl 20>    work    = f->$varName$(iStrt:iLast,{latS:latN},:)   ; (time,lat,lon)<br>ncl 21>    work!0  = "time"                                    ; CAM model names<br>ncl 22>    work!1  = "lat"<br>ncl 23>    work!2  = "lon"<br>ncl 24>    return(work)<br>ncl 25> end<br>ncl 26> ; =========================>  MAIN  <==============================<br>ncl 27> begin<br>ncl 28>    neof    =  2<br>ncl 29> <br>ncl 30>    latS    = -15<br>ncl 31>    latN    =  15<br>ncl 32> <br> ncl 33>    ymdStrt = 20180101                         ; start yyyymmdd<br>ncl 34>    ymdLast = 20181231                         ; last  <br>ncl 35> <br>ncl 36>    yrStrt  = ymdStrt/10000<br>ncl 37>    yrLast  = ymdLast/10000<br>ncl 38> <br>ncl 39>    pltDir  = "/home/peni/"                             ; plot directory<br>*ncl 40>    pltType = "png" <br>ncl 41>    pltName = "mjoclivar"   <br>ncl 42> <br>,ncl 43>    diri    = "/home/peni/"                             ; input directory<br>ncl 44> <br>ncl 45>    filolr  = "<a href="http://olr.day.anomalies.2018.nc" rel="noreferrer noreferrer" target="_blank">olr.day.anomalies.2018.nc</a>"<br>ncl 46>    filu200 = "<a href="http://uwnd.anomalies.nc" rel="noreferrer noreferrer" target="_blank">uwnd.anomalies.nc</a>"<br>ncl 47>    filu850 = "<a href="http://uwnd2.anomalies.nc.nc" rel="noreferrer noreferrer" target="_blank">uwnd2.anomalies.nc.nc</a>"<br>ncl 48> <br>ncl 49> ;************************************************<br>ncl 50> ; create BandPass Filter<br>ncl 51> ;************************************************<br>ncl 52>   ihp      = 2                             ; bpf=>band pass filter<br>ncl 53>   nWgt     = 201<br>ncl 54>   sigma    = 1.0                           ; Lanczos sigma<br>ncl 55>   fca      = 1./100.<br>ncl 56>   fcb      = 1./20.<br>ncl 57>   wgt      = filwgts_lanczos (nWgt, ihp, fca, fcb, sigma )<br>ncl 58> <br>)ncl 59> ;***********************************************************<br>ncl 60> ; Find the indices corresponding to the start/end times<br>ncl 61> ;***********************************************************<br>ncl 62>    f       = addfile (diri+filolr , "r")                         <br>ncl 63>    TIME    = f->time                          ; days since ...<br>ncl 64> <br> ncl 65>    YMD     = cd_calendar(TIME, -2)            ; entire (time,6)<br>ncl 66> <br>ncl 67>    iStrt   = ind(YMD.eq.ymdStrt)              ; index start<br>ncl 68>    iLast   = ind(YMD.eq.ymdLast)              ; index last <br>ncl 69>    delete([/ TIME, YMD /])<br>ncl 70> <br>ncl 71> ;***********************************************************<br>ncl 72> ; Read anomalies<br>ncl 73> ;***********************************************************<br>ncl 74> <br>ncl 75>    work    = read_rename(f,"OLR_anom",iStrt,iLast,latS,latN) ; (time,lat,lon) <br>ncl 76>    OLR     = dim_avg_n_Wrap(work, 1)                         ; (time,lon)<br>ncl 77>    delete(work)<br>ncl 78>    <br>ncl 79>    f       = addfile (diri+filu850 , "r")                         <br>ncl 80>    work    = read_rename(f,"U_anom",iStrt,iLast,latS,latN) ; (time,lat,lon) <br>ncl 81>    U850    = dim_avg_n_Wrap(work, 1)          ; (time,lon)<br>ncl 82>    delete(work)<br>ncl 83>    <br>ncl 84>    f       = addfile (diri+filu200 , "r")                         <br>ncl 85>    work    = read_rename(f,"U_anom",iStrt,iLast,latS,latN) ; (time,lat,lon) <br>ncl 86>    U200    = dim_avg_n_Wrap(work, 1)          ; (time,lon)<br>ncl 87> <br>ncl 88>    dimw    = dimsizes( work )<br>ncl 89>    ntim    = dimw(0)<br>ncl 90>    nlat    = dimw(1)<br>ncl 91>    mlon    = dimw(2)<br>ncl 92>    delete(work)<br>ncl 93> <br>ncl 94>    lon     = OLR&lon                                          <br>ncl 95>    time    = OLR&time                         <br>ncl 96>    date    = cd_calendar(time, -2)            ; yyyymmdd<br>ncl 97> <br>ncl 98> ;************************************************<br>ncl 99> ; Apply the band pass filter to the original anomalies<br>ncl 100> ;************************************************<br>ncl 101>    olr   = wgt_runave_n_Wrap ( OLR, wgt, 0, 0) ; (time,lon)<br>) <br>;************************************************<br>                   <br>  imax_olr_eof1   = maxind(ceof(0,0,:))         <br>  imax_olr_eof2 ncl 102>    u850  = wgt_runave_n_Wrap (U850, wgt, 0, 0)<br>ncl 103>    u200  = wgt_runave_n_Wrap (U200, wgt, 0, 0)<br>ncl 104> <br>ncl 105> ;************************************************<br>ncl 106> ; remove temporal means of band pass series: *not* necessary <br>ncl 107> ;************************************************<br>ncl 108>    olr   = dim_rmvmean_n( olr, 0)              ; (time,lon)<br>ncl 109>    u850  = dim_rmvmean_n(u850, 0)<br>ncl 110>    u200  = dim_rmvmean_n(u200, 0)<br>ncl 111> <br>ncl 112> ;************************************************<br>ncl 113> ; Compute the temporal variance at each lon<br>ncl 114> ;************************************************<br>ncl 115>    var_olr  = dim_variance_n_Wrap( olr, 0)     ; (lon)<br>ncl 116>    var_u850 = dim_variance_n_Wrap(u850, 0)<br>ncl 117>    var_u200 = dim_variance_n_Wrap(u200, 0)<br>ncl 118> <br>ncl 119> ;************************************************<br>ncl 120> ; Compute the zonal mean of the temporal variance<br> Oncl 121> ;************************************************<br>ncl 122>   zavg_var_olr  = dim_avg_n_Wrap( var_olr , 0)     <br>ncl 123>   zavg_var_u850 = dim_avg_n_Wrap( var_u850, 0)<br>ncl 124>   zavg_var_u200 = dim_avg_n_Wrap( var_u200, 0)<br>ncl 125> <br>ncl 126> ;************************************************<br>ncl 127> ; Normalize by sqrt(avg_var*)<br>ncl 128> ;************************************************<br>ncl 129>   olr   =  olr/sqrt(zavg_var_olr )          ; (time,lon)<br>ncl 130>   u850  = u850/sqrt(zavg_var_u850)<br>ncl 131>   u200  = u200/sqrt(zavg_var_u200)<br>ncl 132> <br>cncl 133> ;************************************************<br>ncl 134> ; Combine the normalized data into one variable<br>ncl 135> ;************************************************<br>ncl 136>   cdata     = new ( (/3*mlon,ntim/), typeof(olr), getFillValue(olr))<br>ncl 137>   do ml=0,mlon-1<br>ncl 138>      cdata(ml       ,:) = (/  olr(:,ml) /)<br>ncl 139>      cdata(ml+  mlon,:) = (/ u850(:,ml) /)<br>ncl 140>      cdata(ml+2*mlon,:) = (/ u200(:,ml) /)<br>ncl 141>   end do<br>ncl 142> <br>=ncl 143> ;************************************************<br>ncl 144> ; Compute **combined** EOF; Sign of EOF is arbitrary<br>ncl 145> ;************************************************<br>ncl 146>   eof_cdata    = eofunc(cdata   , neof, False)      ; (neof,3*mlon)<br>ncl 147>   print("==============")<br>ncl 148>   printVarSummary(eof_cdata)<br>ncl 149>   printMinMax(eof_cdata, True)<br>ncl 150> <br>ncl 151>   eof_ts_cdata = eofunc_ts(cdata,eof_cdata,False)   ; (neof,3*ntim)<br>ncl 152>   print("==============")                                  <br>ncl 153>   printVarSummary(eof_ts_cdata)<br>ncl 154>   printMinMax(eof_ts_cdata, True)<br>ncl 155> <br>incl 156> ;************************************************<br>ncl 157> ; For clarity, explicitly extract each variable. Create time series <br>ncl 158> ;************************************************<br>ncl 159> <br>tncl 160>   nvar = 3  ; "olr", "u850", "u200"<br>ncl 161>   ceof = new( (/nvar,neof,mlon/), typeof(cdata), getFillValue(cdata))<br>ncl 162> <br>dncl 163>   do n=0,neof-1<br>ncl 164>      ceof(0,n,:) = eof_cdata(n,0:mlon-1)      ; olr<br>ncl 165>      ceof(1,n,:) = eof_cdata(n,mlon:2*mlon-1) ; u850<br>ncl 166>      ceof(2,n,:) = eof_cdata(n,2*mlon:)       ; u200<br>ncl 167>   end do<br>ncl 168> <br>(ncl 169>   ceof!0   = "var"<br>ncl 170>   ceof!1   = "eof"<br>ncl 171>   ceof!2   = "lon"   <br>ncl 172>   ceof&lon =  olr&lon<br>ncl 173> <br>incl 174>   ceof_ts        = new( (/nvar,neof,ntim/), typeof(cdata), getFillValue(cdata))<br>ncl 175>   ceof_ts(0,:,:) = eofunc_ts_Wrap( olr(lon|:,time|:),ceof(0,:,:),False)   ; (0,neof,ntim)<br>ncl 176>   ceof_ts(1,:,:) = eofunc_ts_Wrap(u850(lon|:,time|:),ceof(1,:,:),False)   ; (1,neof,ntim)<br>ncl 177>   ceof_ts(2,:,:) = eofunc_ts_Wrap(u200(lon|:,time|:),ceof(2,:,:),False)   ; (2,neof,ntim)<br>ncl 178> <br>ncl 179> ;**********************************************t*<br>ncl 180> ; Add code contributed by Marcus N. Morgan, Florida Institute of Technology; Feb 2015<br>ncl 181> ; Calculate % variance (pcv_ )accounted for by OLR, U850 and U200<br>ncl 182> ;************************************************<br>ncl 183> <br>ncl 184>     pcv_eof_olr  = new(neof,typeof(ceof))<br>ncl 185>     pcv_eof_u850 = new(neof,typeof(ceof))<br>ncl 186>     pcv_eof_u200 = new(neof,typeof(ceof))<br>ncl 187>       <br>ncl 188>     do n=0,neof-1<br>ncl 189>        pcv_eof_olr(n)  = avg((ceof(0,n,:)*sqrt(ceof@eval(n)))^2)*100<br>ncl 190>        pcv_eof_u850(n) = avg((ceof(1,n,:)*sqrt(ceof@eval(n)))^2)*100<br>ncl 191>        pcv_eof_u200(n) = avg((ceof(2,n,:)*sqrt(ceof@eval(n)))^2)*100<br>ncl 192>      ;;print("pcv: neof="+(n+1)+":  "+pcv_eof_olr(n)+"  "+pcv_eof_u850(n)+"  "+pcv_eof_u200(n))<br>ncl 193>     end do<br>ncl 194> <br> ncl 195> ;************************************************<br>ncl 196> ; Change sign of EOFs for spatial structures of PC1 and PC2 <br>ncl 197> ; to represent convection over the tropical Indian Ocean and the tropical western Pacific Ocean, respectively <br>ncl 198> ; (Ad hoc approach) <br>ncl 199> ;************************************************<br>ncl 200>                    <br>ncl 201>   imax_olr_eof1   = maxind(ceof(0,0,:))         <br>ncl 202>   imax_olr_eof2   = maxind(ceof(0,1,:))    <br>ncl 203> <br>)ncl 204>   lonmax_eof1 = ceof&lon(imax_olr_eof1)      ; longitude of max value (i.e. suppressed convection)<br>ncl 205>   lonmax_eof2 = ceof&lon(imax_olr_eof2)<br>ncl 206> <br>sncl 207>   if (lonmax_eof1.ge.60 .and. lonmax_eof1.lt.180) then  ; Change the sign of EOF1 <br>ncl 208>       ceof(:,0,:)       = -ceof(:,0,:)                  ; if OLR is positive<br>ncl 209>       ceof_ts(:,0,:)    = -ceof_ts(:,0,:)               ;  over the tropical Indian Ocean<br>ncl 210>       eof_cdata(0,:)    = -eof_cdata(0,:)<br>ncl 211>       eof_ts_cdata(0,:) = -eof_ts_cdata(0,:)<br>ncl 212>   end if<br>ncl 213> <br>encl 214>   if (lonmax_eof2.ge.120 .and. lonmax_eof2.lt.180) then  ; Change the sign of EOF2<br>ncl 215>       ceof(:,1,:)       = -ceof(:,1,:)                   ; if OLR is positive <br>ncl 216>       ceof_ts(:,1,:)    = -ceof_ts(:,1,:)                ; over the tropical western Pacific Ocean<br>ncl 217>       eof_cdata(1,:)    = -eof_cdata(1,:)<br>ncl 218>       eof_ts_cdata(1,:) = -eof_ts_cdata(1,:)<br>ncl 219>   end if<br>ncl 220> <br>ncl 221>   print("==============")<br>rncl 222>   printVarSummary(eof_cdata)<br>ncl 223>   printMinMax(eof_cdata, True)<br>ncl 224> <br>ncl 225> ;************************************************<br>ncl 226> ; Compute cross correlation of each variable's EOF time series at zero-lag<br>ncl 227> ;************************************************<br>ncl 228>   r_olr_u850  = escorc(ceof_ts(0,:,:) , ceof_ts(1,:,:) )  ; (neof)<br>ncl 229>   r_olr_u200  = escorc(ceof_ts(0,:,:) , ceof_ts(2,:,:) )<br>ncl 230>   r_u850_u200 = escorc(ceof_ts(1,:,:) , ceof_ts(2,:,:) )<br>ncl 231> <br>sncl 232>   print("==============")<br>ncl 233>   do n=0,neof-1<br>  ncl 234>      print("neof="+n \<br>ncl 234>           +"  r_olr_u850=" +sprintf("%4.3f",r_olr_u850(n))  \<br>ncl 234>           +"  r_olr_u200=" +sprintf("%4.3f",r_olr_u200(n))  \<br>ncl 234>           +"  r_u850_u200="+sprintf("%4.3f",r_u850_u200(n)) )<br>ncl 235>   end do<br>ncl 236>   print("==============")<br>ncl 237> <br>ncl 238> ;************************************************<br>5N: "+yrStrt+"-"+yrLast <br>  <br>  do ncl 239> ; Compute cross correlation of the multivariate EOF; EOF 1 vs EOF 2<br>ncl 240> ;************************************************<br>ncl 241> <br> ncl 242>   mxlag     = 25<br>ncl 243>   rlag_01   = esccr(eof_ts_cdata(0,:),eof_ts_cdata(1,:), mxlag)   ; (N,mxlag+1)<br>ncl 244>   rlag_10   = esccr(eof_ts_cdata(1,:),eof_ts_cdata(0,:), mxlag)   ; (N,mxlag+1)<br>ncl 245>   ccr_12    = new ( (/2*mxlag+1/), float)    <br>ncl 246> <br>cncl 247>   ccr_12(mxlag:)    = rlag_10(0:mxlag)   <br>ncl 248>   ccr_12(0:mxlag)   = rlag_01(::-1)       ; reverse order<br>ncl 249> ;;print(ccr_12)<br>ncl 250> <br>sncl 251> <br>Pncl 252> ;************************************************<br>ncl 253> ; Normalize the multivariate EOF 1&2 component time series<br>ncl 254> ; Compute (PC1^2+PC2^2): values > 1 indicate "strong" periods<br>ncl 255> ;************************************************<br>ncl 256>   eof_ts_cdata(0,:) = eof_ts_cdata(0,:)/stddev(eof_ts_cdata(0,:))<br>ncl 257>   eof_ts_cdata(1,:) = eof_ts_cdata(1,:)/stddev(eof_ts_cdata(1,:))<br>ncl 258> <br>ncl 259>   mjo_ts_index      = eof_ts_cdata(0,:)^2 + eof_ts_cdata(1,:)^2 <br>ncl 260>   mjo_ts_index_smt  = runave(mjo_ts_index, 91, 0) ; 91-day running mean<br>ncl 261> <br> ncl 262>   nGood   = num(.not.ismissing(mjo_ts_index))     ; # non-missing<br>ncl 263>   nStrong = num(mjo_ts_index .ge. 1.0)<br>ncl 264>   print("nGood="+nGood+"   nStrong="+nStrong+"   nOther="+(nGood-nStrong))<br>ncl 265> <br>rncl 266> ;************************************************<br>ncl 267> ; Write PC results to netCDF for use in another example.<br>ncl 268> ;************************************************<br>ncl 269>   mjo_ts_index!0    = "time"<br>ncl 270>   mjo_ts_index&time = time <br>ncl 271>   mjo_ts_index@long_name = "MJO PC INDEX" <br>ncl 272>   mjo_ts_index@info      = "(PC1^2 + PC2^2)" <br>ncl 273> <br> ncl 274>   PC1           = eof_ts_cdata(0,:)<br>ncl 275>   PC1!0         = "time"<br>ncl 276>   PC1&time      =  time<br>ncl 277>   PC1@long_name = "PC1"<br>ncl 278>   PC1@info      = "PC1/stddev(PC1)"<br>ncl 279> <br>oncl 280>   PC2           = eof_ts_cdata(1,:)<br>ncl 281>   PC2!0         = "time"<br>ncl 282>   PC2&time      =  time<br>ncl 283>   PC2@long_name = "PC2"<br>ncl 284>   PC2@info      = "PC2/stddev(PC2)"<br>ncl 285> <br>ncl 286>   diro = "./"<br>ncl 287>   filo = "MJO_PC_INDEX.nc"<br>ncl 288>   system("/bin/rm -f "+diro+filo)   ; remove any pre-existing file<br>ncl 289>   ncdf = addfile(diro+filo,"c")     ; open output netCDF file<br>of legend.<br>  ncl 290>                                     ; make time an UNLIMITED dimension <br>ncl 291>   filedimdef(ncdf,"time",-1,True)   ; recommended  for most applications<br>ncl 292>                                     ; output variables directly<br>ncl 293>   ncdf->MJO_INDEX = mjo_ts_index    <br>ncl 294>   ncdf->PC1       = PC1     <br>ncl 295>   ncdf->PC2       = PC2     <br>ncl 296> <br>ncl 297> ;------------------------------------------------------------<br>ncl 298> ; PLOTS<br>ncl 299> ;------------------------------------------------------------<br>ncl 300> <br>ncl 301>   yyyymmdd = cd_calendar(time, -2)<br>tncl 302>   yrfrac   = yyyymmdd_to_yyyyfrac(yyyymmdd, 0.0)<br>ncl 303>   delete([/ yrfrac@long_name, lon@long_name /])<br>ncl 304> <br>ncl 305>   day      = ispan(-mxlag, mxlag, 1)<br>ncl 306>  ;day@long_name = "lag (day)"<br>ncl 307> <br>ncl 308>   pltPath = pltDir+pltName<br>ncl 309> <br>ncl 310>   wks = gsn_open_wks(pltType,pltPath)<br>ncl 311>   gsn_define_colormap(wks,"default")<br>ncl 312>   plot = new(3,graphic)                <br>ncl 313> <br>ncl 314> ;************************************************<br>ncl 315> ; Multivariate EOF plots<br>ncl 316> ;************************************************<br>ncl 317>   rts           = True<br>ncl 318>   rts@gsnDraw   = False       ; don't draw yet<br>ncl 319>   rts@gsnFrame  = False       ; don't advance frame yet<br>ncl 320>   rts@gsnScale  = True        ; force text scaling               <br>ncl 321> <br>ncl 322>   rts@vpHeightF = 0.40        ; Changes the aspect ratio<br>ncl 323>   rts@vpWidthF  = 0.85<br>ncl 324>   rts@vpXF      = 0.10        ; change start locations<br>ncl 325>   rts@vpYF      = 0.75        ; the plot<br>ncl 326>   rts@xyLineThicknesses = (/2, 2, 2/)<br>ncl 327>   rts@xyLineColors      = (/"black","red","blue"/)<br>ncl 328>   rts@gsnYRefLine       = 0.                  ; reference line   <br>ncl 329>   rts@trXMaxF           = max(lon)<br>ncl 330>   rts@trXMinF           = min(lon)<br>ncl 331> <br>ncl 332>   rts@pmLegendDisplayMode    = "Always"            ; turn on legend<br>ncl 333>   rts@pmLegendSide           = "Top"               ; Change location of <br>ncl 334>   rts@pmLegendParallelPosF   = 1.16                ; move units right<br>ncl 335>   rts@pmLegendOrthogonalPosF = -0.50               ; move units down<br>ncl 336>   rts@pmLegendWidthF         = 0.15                ; Change width and<br>ncl 337>   rts@pmLegendHeightF        = 0.15                ; height of legend.<br>ncl 338>   rts@lgLabelFontHeightF     = 0.0175<br>ncl 339> <br>ncl 340> <br>ncl 341>   rtsP                       = True                ; modify the panel plot<br>ncl 342> ;  rtsP@gsnMaximize           = True                ; large format<br>ncl 343>   rtsP@gsnPanelMainString     = "Multivariate EOF: 15S-15N: "+yrStrt+"-"+yrLast <br>ncl 344>   <br>ncl 345>   do n=0,neof-1<br>ncl 346>     rts@xyExplicitLegendLabels = (/"OLR: "+sprintf("%4.1f", pcv_eof_u200(n)) +"%" \<br>ncl 346>                                   ,"U850: "+sprintf("%4.1f", pcv_eof_u850(n))+"%" \<br>ncl 346>                                   ,"U200: "+sprintf("%4.1f", pcv_eof_olr(n))+"%" /)<br>ncl 347>     rts@gsnLeftString  = "EOF "+(n+1)<br>ncl 348>     rts@gsnRightString = sprintf("%3.1f",ceof@pcvar(n))  +"%"<br>ncl 349>     plot(n) = gsn_csm_xy (wks,lon,ceof(:,n,:),rts)<br>ncl 350>   end do<br>ncl 351>   gsn_panel(wks,plot(0:1),(/2,1/),rtsP)     ; now draw as one plot<br>ncl 352> <br>ncl 353> ;-----------------------------------------<br>ncl 354> ; The following doesn't work with some older versions of NCL<br>ncl 355> ; With old versions, the user must delete each individually.<br>ncl 356> ;----------------------------------------- <br>ncl 357>   delete([/ rts@xyExplicitLegendLabels, rts@pmLegendDisplayMode \<br>ncl 357>           , rts@xyLineThicknesses     , rts@gsnLeftString       \<br>ncl 357>           , rts@gsnRightString        , rts@xyLineColors        \<br>ncl 357>           , rts@trXMaxF               , rts@trXMinF             /] )<br>ncl 358> <br>ncl 359>   lag                        = ispan(-mxlag,mxlag,1)<br>ncl 360>   lag@long_name              = "lag (days)"<br>ncl 361> <br>ncl 362>   plot(0)                    = gsn_csm_xy (wks, lag ,ccr_12,rts)<br>ncl 363>   rtsP@gsnPanelMainString    = "Cross Correlation: Multivariate EOF: 15S-15N: " \<br>ncl 363>                    +  yrStrt+"-"+yrLast <br>ncl 364>   rtsP@gsnPaperOrientation   = "portrait"        ; force portrait<br>ncl 365>   gsn_panel(wks,plot(0),(/1,1/),rtsP)     ; now draw as one plot<br>ncl 366> <br>ncl 367> ;************************************************<br>ncl 368> ; MJO "strong" index <br>ncl 369> ;************************************************<br>ncl 370>   rts@gsnYRefLine        = 1.0<br>ncl 371>   rts@gsnYRefLineColor   = "black"<br>ncl 372>   rts@xyMonoDashPattern  = True<br>ncl 373>   rts@xyLineColors       = (/"black", "blue"/)<br>ncl 374>   rts@xyLineThicknesses  = (/1, 2/)<br>ncl 375>   rts@pmLegendDisplayMode    = "Always"            ; turn on legend<br>ncl 376>   rts@pmLegendWidthF         = 0.12                ; Change width and<br>ncl 377>   rts@pmLegendHeightF        = 0.10                ; height of legend.<br>ncl 378>   rts@pmLegendParallelPosF   = 0.86                ; move units right<br>ncl 379>   rts@pmLegendOrthogonalPosF = -0.40               ; move units down<br>ncl 380>   rts@xyExplicitLegendLabels = (/"daily", "91-day runavg" /)<br>ncl 381> <br>ncl 382>   mjo_ind_plt = new ( (/2,ntim/), typeof(mjo_ts_index))<br>ncl 383>   mjo_ind_plt(0,:) = mjo_ts_index<br>ncl 384>   mjo_ind_plt(1,:) = (/ mjo_ts_index_smt /)<br>ncl 385>   plot(0) = gsn_csm_xy(wks, yrfrac,mjo_ind_plt,rts)<br>ncl 386> <br>ncl 387>   rtsP@gsnPanelMainString   = "MJO Index: (PC1^2+ PC2^2) : 15S-15N: "+yrStrt+"-"+yrLast <br>ncl 388>   gsn_panel(wks,plot(0),(/1,1/),rtsP)     ; now draw as one plot<br>ncl 389> <br>ncl 390>  end<br><div><br></div><div>Thank you<br></div><br></div>
_______________________________________________<br>
ncl-talk mailing list<br>
<a href="mailto:ncl-talk@ucar.edu" rel="noreferrer noreferrer" target="_blank">ncl-talk@ucar.edu</a><br>
List instructions, subscriber options, unsubscribe:<br>
<a href="http://mailman.ucar.edu/mailman/listinfo/ncl-talk" rel="noreferrer noreferrer noreferrer" target="_blank">http://mailman.ucar.edu/mailman/listinfo/ncl-talk</a></blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div><div><div><span><font color="#888888">Adam Phillips <br></font></span></div><span><font color="#888888">Associate Scientist,  </font></span><span><font color="#888888">Climate and Global Dynamics Laboratory, NCAR<br></font></span></div></div><div><span><font color="#888888"><a href="http://www.cgd.ucar.edu/staff/asphilli/" rel="noreferrer noreferrer" target="_blank">www.cgd.ucar.edu/staff/asphilli/</a>   </font></span><span><font color="#888888">303-497-1726 </font></span></div><span><font color="#888888"></font></span><div><div><span><font color="#888888"><br></font></span><div><span><font color="#888888"><a href="http://www.cgd.ucar.edu/staff/asphilli" rel="noreferrer noreferrer" target="_blank"></a></font></span></div></div></div></div></div></div></div></div></div></div></div>
</blockquote></div></div>
_______________________________________________<br>
ncl-talk mailing list<br>
<a href="mailto:ncl-talk@ucar.edu" rel="noreferrer noreferrer" target="_blank">ncl-talk@ucar.edu</a><br>
List instructions, subscriber options, unsubscribe:<br>
<a href="http://mailman.ucar.edu/mailman/listinfo/ncl-talk" rel="noreferrer noreferrer noreferrer" target="_blank">http://mailman.ucar.edu/mailman/listinfo/ncl-talk</a></blockquote></div>
</blockquote></div>
</blockquote></div>
</div>
</div>
</blockquote></div>
_______________________________________________<br>
ncl-talk mailing list<br>
<a href="mailto:ncl-talk@ucar.edu" target="_blank">ncl-talk@ucar.edu</a><br>
List instructions, subscriber options, unsubscribe:<br>
<a href="http://mailman.ucar.edu/mailman/listinfo/ncl-talk" rel="noreferrer" target="_blank">http://mailman.ucar.edu/mailman/listinfo/ncl-talk</a></blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr"><div dir="ltr"><div><div dir="ltr">Barry H. Lynn, Ph.D<div><div>Senior Associate Scientist, Lecturer,</div><div><div><span style="color:rgb(136,136,136)">The Institute of the Earth Science, </span><br style="color:rgb(136,136,136)"><span style="color:rgb(136,136,136)">The Hebrew University of Jerusalem, </span><br style="color:rgb(136,136,136)"><span style="color:rgb(136,136,136)">Givat Ram, Jerusalem 91904, Israel </span><br style="color:rgb(136,136,136)"></div><span style="color:rgb(136,136,136)">Tel: 972 547 231 170</span><br style="color:rgb(136,136,136)"><span style="color:rgb(136,136,136)">Fax: (972)-25662581</span></div></div><div><span style="color:rgb(136,136,136)"><br></span></div><div>C.E.O, Weather It Is, LTD<br>Weather and Climate Focus<br><a href="http://weather-it-is.com" target="_blank">http://weather-it-is.com</a><br>Jerusalem, Israel<br>Local: 02 930 9525<br>Cell: 054 7 231 170<br>Int-IS: x972 2 930 9525<br><br></div></div></div></div></div>
</blockquote></div>
_______________________________________________<br>
ncl-talk mailing list<br>
<a href="mailto:ncl-talk@ucar.edu" target="_blank">ncl-talk@ucar.edu</a><br>
List instructions, subscriber options, unsubscribe:<br>
<a href="http://mailman.ucar.edu/mailman/listinfo/ncl-talk" rel="noreferrer" target="_blank">http://mailman.ucar.edu/mailman/listinfo/ncl-talk</a></blockquote></div>