[ncl-talk] Sub: error

Adv advita6 at gmail.com
Mon Jun 6 16:54:03 MDT 2016


Hi ,
Could someone help me to fix this bug? I don't whats the fatal error says
about.
Thank you

;***************************************************
begin
;***************************************************
  yrStrt = 196501
  yrLast = 200512

  season = "SON"    ; choose Dec-Jan-Feb seasonal mean
;####################################################
............
f    = addfile("hgt.mon.mean.nc", "r")   ; note the "s" of addfile
u = f->hgt
printVarSummary (u)
   YYYYMM = cd_calendar( f->time, -1)

   iStrt  = ind(YYYYMM.eq.yrStrt)
   iLast  = ind(YYYYMM.eq.yrLast)
n=u(level|:,lat|:,lon|:,time|iStrt:iLast)
anew=n({500},:,:,:)
printVarSummary (anew)
;############################################################
;anomalies-------------------------
temp= rmMonAnnCycTLL(anew(time|:,lat|:,lon|:))
printVarSummary(temp)
;dtrend the values/------------------
anew1Dtrend = dtrend(temp(lat|:,lon|:,time|:),True)
printVarSummary(anew1Dtrend)
;print(anew1Dtrend at slope)
;print(anew1Dtrend at y_intercept)
copy_VarCoords(anew,anew1Dtrend)

; ==============================================================
; compute desired global seasonal mean: month_to_season (contributed.ncl)
; ==============================================================
  anew1a    = month_to_season (anew1Dtrend(time|:,lat|:,lon|:), season)
  nyrs   = dimsizes(anew1a&time)
  printVarSummary(anew1a)
printMinMax(anew1a,False)
;%%%%%%%%%%%%%%%EOF%%%%%%%%%%%%%%%%%%%%%%%%
;Dont change at all . Here is a trick..
;Globe
  latS   = 37
  latN   =  49.
 lonL     =  -116
 lonR     = -90

;Center=0

;%%%%%%%%%%%%%%%EOF%%%%%%%%%%%%%%%%%%%%%%%%


neof   = 3        ; number of EOFs
  optEOF = True
  optEOF at jopt = 0   ; This is the default; most commonly used; no need to
specify.
;;optEOF at jopt = 1   ; **only** if the correlation EOF is desired
; ==============================================================
; dataset longitudes span 0=>357.5
; Because EOFs of the North Atlantic Oscillation are desired
; use the "lonFlip" (contributed.ncl) to reorder
; longitudes to span -180 to 177.5: facilitate coordinate subscripting
; ==============================================================
  anew1a   = lonFlip(anew1a)
  printVarSummary(anew1a)                              ; note the longitude
coord
orig=anew1a({lat|latS:latN},{lon|lonL:lonR},time|:)
  printVarSummary(orig)

; =================================================================
; create weights:  sqrt(cos(lat))   [or sqrt(gw) ]
; =================================================================
  rad    = 4.*atan(1.)/180.
  clat   = anew1a(0,:,0)
  printVarSummary(clat)
  clat   = sqrt( cos(rad*clat) )                 ; gw for gaussian grid
  printVarSummary(clat)
; =================================================================
; weight all observations
; =================================================================
  wSLP   = anew1a                                   ; copy meta data
  printVarSummary(wSLP)
  printMinMax(wSLP,False)
  wSLP   = anew1a*conform(anew1a, clat, 1)
  wSLP at long_name = "Wgt: "+anew1a at long_name
  printVarSummary(wSLP)
printMinMax(wSLP,False)
; =================================================================
; Reorder (lat,lon,time) the *weighted* input data
; Access the area of interest via coordinate subscripting
; =================================================================
  x      = wSLP({lat|latS:latN},{lon|lonL:lonR},time|:)
  printVarSummary(x)
print(min(x))
print(max(x))
 dimx   = dimsizes( x )
  mln    = dimx(1)
  sumWgt = mln*sum( clat({lat|latS:latN}) )
printMinMax(sumWgt,False)
  eof    = eofunc_Wrap(x, neof, optEOF)
  printVarSummary(eof)

print(eof(0,:,:))
;return
;;;;EOF spatial patterns may be tested for orthogonality by using the dot
product:
d01 = sum(eof(0,:,:)*eof(1,:,:))
  d12 = sum(eof(1,:,:)*eof(2,:,:))
  d02 = sum(eof(0,:,:)*eof(2,:,:))
  print("d01="+d01+"  d12="+d12+"  d02="+d02)  ; may be +/- 1e-8
;return
printMinMax(eof(0,:,:),False)
;To retain same unit
  do ne=0,neof-1
     eof(ne,:,:) = eof(ne,:,:)*sqrt( eof at eval(ne) )
  end do
;print(eof(:,:,:))
printMinMax(eof(0,:,:),False)
;return
;========================================================
;Print the Percentage variance;;;;;;
; print("% var="+ eof at pcvar(0))
;========================================================
eof_ts = eofunc_ts_Wrap (x, eof, False)
 printVarSummary( eof_ts )
  print( eof_ts )
;=======================================================
;To test the EOF time series for orthogonality, compute correlations;;;;;
  r01 = escorc(eof_ts(0,:), eof_ts(1,:))
  r12 = escorc(eof_ts(1,:), eof_ts(2,:))
  r02 = escorc(eof_ts(0,:), eof_ts(2,:))
  print("r01="+r01+"  r12="+r12+"  r02="+r02)   ; numbers may be +/- 1e-8
; =================================================================
; Normalize time series: Sum spatial weights over the area of used
; =================================================================
;  dimx   = dimsizes( x )
;  mln    = dimx(1)
;  sumWgt = mln*sum( clat({lat|latS:latN}) )
;eof_ts = eof_ts/sumWgt
;nlat = dimsizes(eof&lat)    ; lat==> whatever the coordinate name is
;  mlon= dimsizes(eof&lon)
;print(eof_ts)
  eof_ts = dim_standardize_n( eof_ts, 0, 1)      ; normalize
  eof_regres = eof                               ; create an array w meta
data
  do ne=0,neof-1
     eof_regres(ne,:,:) = (/ regCoef(eof_ts(ne,:),orig(lat|:,lon|:,time|:))
/)
  end do
printMinMax(eof_regres(0,:,:),False)
printVarSummary( eof_regres )
;print(eof_ts)
;print(eof(0,:,:))
;===================================================================
  nlat = dimsizes(eof&lat)    ; lat==> whatever the coordinate name is
  mlon= dimsizes(eof&lon)
;------------------------------------------------------------

; =================================================================
; Extract the YYYYMM from the time coordinate
; associated with eof_ts [same as SLP&time]
; =================================================================

  yyyymm = cd_calendar(eof_ts&time,-2)/100
;;yrfrac = yyyymm_to_yyyyfrac(yyyymm, 0.0); not used here
;*******************************************
; North significance test: any pcvar, eval_transpose, eval can be used
;*******************************************

  dimp   = dimsizes(x)
  ntim   = dimp(0)                                            ; max #
eigenvalues possible

  prinfo = True
  sig    = eofunc_north(eof at pcvar, ntim, prinfo)

;************************************************
; plotting parameters
;************************************************
;wks =
gsn_open_wks("x11","SlpReconst_Anom_eof_Reanal_"+season+"_temp1965-2005")
wks =
gsn_open_wks("x11","heightReconst_Anom_eof_Reanal_"+season+"_temp1965-2005")
 gsn_define_colormap(wks,"NCV_jaisnd")
 plot = new(neof,graphic)                ; create graphic array
                                          ; only needed if paneling
; EOF patterns

  res                      = True
  res at gsnDraw              = False        ; don't draw yet
  res at gsnFrame             = False        ; don't advance frame yet

;---This resource not needed in V6.1.0
  res at gsnSpreadColors      = True         ; spread out color table
res at cnLevelSelectionMode = "ManualLevels"  ; set manual contour levels
  res at cnMinLevelValF       = -10.0       ; set min contour level
  res at cnMaxLevelValF       = 10.0   ; set max contour level
  res at cnLevelSpacingF      = 0.5      ; set contour spacing

  res at gsnAddCyclic         = False        ; plotted dataa are not cyclic

  res at mpFillOn             = False        ; turn off map fill

  res at mpMinLatF            = eof&lat(0)+3         ; zoom in on map
  res at mpMaxLatF            = eof&lat(nlat-1)
  res at mpMinLonF            = eof&lon(0)
  res at mpMaxLonF            = eof&lon(mlon-1)

res at cnLineLabelsOn       = False
 res at cnFillOn             = True         ; turn on color fill
  res at cnLinesOn            = False        ; True is default
  res at lbLabelBarOn         = False        ; turn off individual lb's
; res at mpPerimOn              = True                    ; draw box around map
res at mpGeophysicalLineThicknessF = 3.0
res at mpGeophysicalLineColor = "Black"; (/22/)
res at mpOutlineBoundarySets = "GeophysicalAndUSStates" ; add state boundaries
res at mpNationalLineColor  = res at mpGeophysicalLineColor
res at mpUSStateLineThicknessF = 3.0
res at mpUSStateLineColor  = res at mpGeophysicalLineColor


                                          ; set symmetric plot min/max
;  symMinMaxPlt(eof, 16, False, res)       ; contributed.ncl

; panel plot only resources
  resP                     = True         ; modify the panel plot
  resP at gsnMaximize         = True         ; large format
  resP at gsnPanelLabelBar    = True         ; add common colorbar
  resP at lbLabelAutoStride   = True         ; auto stride on labels
  resP at lbLabelFontHeightF      = 0.017     ; label bar font

  yStrt                    = yyyymm(0)/100
  yLast                    = yyyymm(nyrs-1)/100
  resP at txString            = "Regres_Reanal_height(m): "+season+":
"+yStrt+"-"+yLast

;*******************************************
; first plot
;*******************************************
  do n=0,neof-1
     res at gsnLeftStringFontHeightF   = 0.02        ; instead of using
txFontHeightF or gsnStringFontHeightF
  res at gsnCenterStringFontHeightF = 0.02            ; to set the
gsnLeft/Center/RightString font heights,
  res at gsnRightStringFontHeightF  = 0.02        ; individually set each
string's font height.
     res at gsnLeftString  = "EOF "+(n+1)
     res at gsnCenterString= "NORTH="+sig(n)
     res at gsnRightString = sprintf("%5.1f", eof at pcvar(n)) +"%"
     plot(n)=gsn_csm_contour_map_ce(wks,eof_regres(n,:,:),res)
  end do
  gsn_panel(wks,plot,(/neof,1/),resP)               ; now draw as one plot


  frame(wks)


end



ERROR printed:

fatal:Loop variable must be scalar, can't execute loop
fatal:["Execute.c":8575]:Execute: Error occurred at or near line 242 in
file heightReconstruct_Temp_eof_Reanal.ncl
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160606/4fe1c40b/attachment.html 


More information about the ncl-talk mailing list