[ncl-talk] Subscript out of range, error in subscript #0
ANDIKA FAUZIAH HAPSARI, A.P
andika.hapsari at bmkg.go.id
Thu Dec 5 20:16:17 MST 2019
Hi NCL expert
I try to filter Kelvin waves in U and V wind data from ECMWF. But I got this
error..
fatal:["NclFile.c":2120]:Subscript out of range, error in subscript #0
fatal:["Execute.c":8637]:Execute: Error occurred at or near line 102 in file
filter_waves4_uwnd.ncl
Does anyone know what the error means and what should I do? Thanks in
advance..
Here is Carl Shreck's script that I have modified
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; filter_waves.ncl
; Carl Schreck (carl at cicsnc.org)
; February 2011
; Updated for style October 2011
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Description: Filter data for each equatorial wave type
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
;load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/diagnostics_cam.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/kf_filter.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/cd_string.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/print_clock.ncl"
;load "$CJS_NCL_LIB/print_clock.ncl"
begin
print_clock( "Here we go! " )
; These are some parameters that could be useful to have up top
mis = -999
mis at _FillValue = -999
; xmin = 70
; xmax = 170
ymin = -15
ymax = 15
zmin = 1000
zmax = 10
units = "hours since 1900-01-01 00:00:00.0"
tmin = cd_inv_calendar( 2018, 12, 01, 00, 0, 0, units, 0 )
tmax = cd_inv_calendar( 2018, 12, 31, 00, 0, 0, units, 0 )
if( .not.isvar("filtName") ) then
filtName = "k09"
end if
obsPerDay = 4
if( .not.isvar("varName") ) then
varName = "u"
end if
print( varName )
basePath = "/home/sarihmmm/datareanly/reanalysis/ecmwf/"
; basePath = "~/data/nasa_ymc/fcst_verif/"
; basePath = "~/data/olr/current/"
; pathIn = basePath + ".std.nc"
; pathOut = basePath + ".waves.std.nc"
pathIn = basePath + "UVwind_dec2018.nc"
pathOut = basePath + "anginU_dec2018.kelvin.nc"
; pathOut = basePath + "trmm3b42." + filtName + ".nc"
; pathIn = basePath + varName + ".anom.nc"
; pathOut = basePath + varName + "." + filtName + ".nc"
; pathIn = basePath + "anom/" + varName + ".anom.nc"
; pathOut = basePath + "waves/" + varName + ".anom.waves.nc"
calcHigh = False
calcLow = False
calcMrg = False
calcMjo = False
calcKelvin = True
calcEr = False
calcMtd = False
calcTd = False
calcWig = False
makeNewFile = True
; Open the input files
fin = addfile( pathIn, "r" )
longitude = fin->longitude
latitude = fin->latitude({ymin:ymax})
; level = fin->level
level = fin->level({zmin:zmax})
time = fin->time({tmin:tmax})
latitude at actual_range = (/ min(latitude), max(latitude) /)
; Open the output files
setfileoption("nc","Format","LargeFile")
; setfileoption("nc","Format","NetCDF4")
; setfileoption("nc","CompressionLevel",1)
if( makeNewFile ) then
system( "rm " + pathOut )
fout = addfile( pathOut, "c" )
fout->time = time
fout->latitude = latitude
fout->longitude = longitude
else
fout = addfile( pathOut, "w" )
end if
print_clock( "Reading the input data... " )
data = tofloat(
fin->$varName$(time({tmin:tmax}),{zmin:zmax},{ymin:ymax},:) )
data&time at beginning_date = cd_string( data&time(0), "" )
data&time at ending_date = cd_string( data&time(dimsizes(data&time)-1), ""
)
if( any( ismissing(data) ) ) then
print_clock( "WARNING: Setting missing data to zero" )
data = where( ismissing(data), 0, data )
end if
printVarSummary(data)
fout->unfilt = data
if( calcHigh ) then
print_clock( "Filtering High..." )
high = data
high at long_name = "High frequency"
high at filter = "Highpass time filter"
high at wavenumber = (/ -9999, 9999 /)
high at period = (/ 0.01, 120 /)
high at depth = (/ mis, mis /)
do y = 0, ( dimsizes(lat) - 1 )
; print_clock( y + " " + ( dimsizes(lat) - 1 ) + " " )
high(time|:,{lat|lat(y)},lon|:) = (/ kf_filter( \\
data(time|:,{lat|lat(y)},lon|:), obsPerDay, \\
high at period(0), high at period(1), \\
high at wavenumber(0), high at wavenumber(1), \\
high at depth(0), high at depth(1), \\
"none" ) /)
end do
print_clock( "Writing high... " )
fout->high = high
delete(high)
end if
if( calcLow ) then
print_clock( "Filtering Low..." )
low = data
low at long_name = "Low frequency"
low at filter = "Lowpass time filter"
if( filtName.eq."wk99" ) then
low at wavenumber = (/ -99, 99 /)
low at period = (/ 120, 9999 /)
low at depth = (/ mis, mis /)
else
low at wavenumber = (/ -10, 10 /)
low at period = (/ 120, 9999 /)
low at depth = (/ mis, mis /)
end if
do y = 0, ( dimsizes(lat) - 1 )
; print_clock( y + " " + ( dimsizes(lat) - 1 ) + " " )
low(time|:,{lat|lat(y)},lon|:) = (/ kf_filter( \\
data(time|:,{lat|lat(y)},lon|:), obsPerDay, \\
low at period(0), low at period(1), \\
low at wavenumber(0), low at wavenumber(1), \\
low at depth(0), low at depth(1), \\
"none" ) /)
end do
print_clock( "Writing Low... " )
fout->low = low
; delete(low)
end if
if( calcMrg ) then
print_clock( "Filtering MRG..." )
mrg = data
mrg at long_name = "Mixed Rossby-Gravity Waves in " +
str_upper(varName)
mrg at filter = "Wheeler & Kiladis (1999)"
mrg at wavenumber = (/ -10, -1 /)
mrg at period = (/ 3, 96 /)
mrg at depth = (/ 8, 90 /)
do y = 0, ( dimsizes(lat) - 1 )
; print_clock( y + " " + ( dimsizes(lat) - 1 ) + " " )
mrg(time|:,{lat|lat(y)},lon|:) = (/ kf_filter( \\
data(time|:,{lat|lat(y)},lon|:), obsPerDay, \\
mrg at period(0), mrg at period(1), \\
mrg at wavenumber(0), mrg at wavenumber(1), \\
mrg at depth(0), mrg at depth(1), \\
"mrg" ) /)
end do
print_clock( "Writing MRG... " )
fout->mrg = mrg
delete(mrg)
end if
if( calcMjo ) then
print_clock( "Filtering MJO..." )
mjo = data
mjo at long_name = "Madden-Julian Oscillation in " + str_upper(varName)
if( filtName.eq."wk99" ) then
mjo at filter = "Wheeler & Kiladis (1999)"
mjo at wavenumber = (/ 1, 5 /)
mjo at period = (/ 30, 96 /)
mjo at depth = (/ mis, mis /)
else
mjo at filter = "Kiladis et al. (2005 JAS) for 20-100"
mjo at wavenumber = (/ 0, 9 /)
mjo at period = (/ 20, 100 /)
mjo at depth = (/ mis, mis /)
end if
do y = 0, ( dimsizes(lat) - 1 )
; print_clock( y + " " + ( dimsizes(lat) - 1 ) + " " )
mjo(time|:,{lat|lat(y)},lon|:) = (/ kf_filter( \\
data(time|:,{lat|lat(y)},lon|:), obsPerDay, \\
mjo at period(0), mjo at period(1), \\
mjo at wavenumber(0), mjo at wavenumber(1), \\
mjo at depth(0), mjo at depth(1), \\
"none" ) /)
end do
print_clock( "Writing MJO... " )
fout->mjo = mjo
sumOfModes = data
sumOfModes = low + mjo
sumOfModes at long_name = "Low+MJO"
fout->low_mjo = sumOfModes
; delete(mjo)
end if
if( calcEr ) then
print_clock( "Filtering ER..." )
er = data
er at long_name = "Equatorial Rossby Waves in " + str_upper(varName)
if( filtName.eq."wk99" ) then
er at filter = "Wheeler & Kiladis (1999)"
er at wavenumber = (/ -10, -1 /)
er at period = (/ 9.7, 48 /)
er at depth = (/ 8, 90 /)
else
er at filter = "Kiladis et al. (2009 Rev. Geophys.)"
er at wavenumber = (/ -10, -1 /)
er at period = (/ 9.7, 72 /)
if( filtName.eq."k09" ) then
er at depth = (/ mis, 90 /)
else ; wide
er at depth = (/ mis, mis /)
end if
end if
do y = 0, ( dimsizes(lat) - 1 )
; print_clock( y + " " + ( dimsizes(lat) - 1 ) + " " )
er(time|:,{lat|lat(y)},lon|:) = (/ kf_filter( \\
data(time|:,{lat|lat(y)},lon|:), obsPerDay, \\
er at period(0), er at period(1), \\
er at wavenumber(0), er at wavenumber(1), \\
er at depth(0), er at depth(1), \\
"er" ) /)
end do
print_clock( "Writing ER... " )
fout->er = er
sumOfModes = sumOfModes + er
sumOfModes at long_name = "Low+MJO+ER"
fout->low_mjo_er = sumOfModes
; delete(er)
end if
if( calcKelvin ) then
print_clock( "Filtering KELVIN..." )
kelvin = data
kelvin at long_name = "Kelvin Waves in " + str_upper(varName)
if( filtName.eq."wk99" ) then
kelvin at filter = "Wheeler & Kiladis (1999)"
kelvin at wavenumber = (/ 1, 14 /)
kelvin at period = (/ 2.5, 30 /)
kelvin at depth = (/ 8, 90 /)
else if( filtName.eq."k09" ) then
kelvin at filter = "Straub & Kiladis (2002) to 20 days"
kelvin at wavenumber = (/ 1, 14 /)
kelvin at period = (/ 2.5, 20 /)
kelvin at depth = (/ 8, 90 /)
else if( filtName.eq."wide" ) then
kelvin at filter = "Broader square filter"
kelvin at wavenumber = (/ 1, 14 /)
kelvin at period = (/ 2.5, 20 /)
kelvin at depth = (/ mis, mis /)
end if
end if
end if
do y = 0, ( dimsizes(latitude) - 1 )
; print_clock( y + " " + ( dimsizes(latitude) - 1 ) + " " )
kelvin(time|:,{latitude|latitude(y)},longitude|:) = (/ kf_filter( \\
data(time|:,{latitude|latitude(y)},longitude|:), obsPerDay, \\
kelvin at period(0), kelvin at period(1), \\
kelvin at wavenumber(0), kelvin at wavenumber(1), \\
kelvin at depth(0), kelvin at depth(1), \\
"kelvin" ) /)
end do
print_clock( "Writing KELVIN... " )
fout->kelvin = kelvin
; sumOfModes = sumOfModes + kelvin
; sumOfModes at long_name = "Low+MJO+ER+Kelvin"
; fout->low_mjo_er_kelvin = sumOfModes
; delete(kelvin)
end if
if( calcMtd ) then
print_clock( "Filtering MRG/TD..." )
td = data
td at long_name = "MRG/TD in " + str_upper(varName)
td at filter = "Frank and Roundy (2006), extended to wavenumber
-20"
td at wavenumber = (/ -20, -0 /)
td at period = (/ 2.5, 10 /)
td at depth = (/ mis, mis /)
do y = 0, ( dimsizes(lat) - 1 )
; print_clock( y + " " + ( dimsizes(lat) - 1 ) + " " )
td(time|:,{lat|lat(y)},lon|:) = (/ kf_filter( \\
data(time|:,{lat|lat(y)},lon|:), obsPerDay, \\
td at period(0), td at period(1), \\
td at wavenumber(0), td at wavenumber(1), \\
td at depth(0), td at depth(1), \\
"none" ) /)
end do
print_clock( "Writing TD... " )
fout->td = td
delete(td)
end if
if( calcTd ) then
print_clock( "Filtering TD..." )
td = data
td at long_name = "TD-type Disturbances in " + str_upper(varName)
td at filter = "Roundy and Frank (2004), extended to wavenumber
-20"
td at wavenumber = (/ -20, -6 /)
td at period = (/ 2.5, 5 /)
td at depth = (/ 90, mis /)
do y = 0, ( dimsizes(lat) - 1 )
; print_clock( y + " " + ( dimsizes(lat) - 1 ) + " " )
td(time|:,{lat|lat(y)},lon|:) = (/ kf_filter( \\
data(time|:,{lat|lat(y)},lon|:), obsPerDay, \\
td at period(0), td at period(1), \\
td at wavenumber(0), td at wavenumber(1), \\
td at depth(0), td at depth(1), \\
"none" ) /)
end do
print_clock( "Writing TD... " )
fout->td = td
delete(td)
end if
if( calcWig ) then
print_clock( "Filtering WIG..." )
wig = data
wig at long_name = "Westward Inertial-Gravity Waves in " +
str_upper(varName)
wig at filter = "Wheeler & Kiladis (1999), n=1"
wig at wavenumber = (/ -15, -1 /)
wig at period = (/ 1.25, 3 /)
wig at depth = (/ 12, 90 /)
do y = 0, ( dimsizes(lat) - 1 )
; print_clock( y + " " + ( dimsizes(lat) - 1 ) + " " )
wig(time|:,{lat|lat(y)},lon|:) = (/ kf_filter( \\
data(time|:,{lat|lat(y)},lon|:), obsPerDay, \\
wig at period(0), wig at period(1), \\
wig at wavenumber(0), wig at wavenumber(1), \\
wig at depth(0), wig at depth(1), \\
"ig1" ) /)
end do
print_clock( "Writing WIG... " )
fout->wig = wig
delete(wig)
end if
print_clock( "Closing file..." )
delete(fout)
print_clock( "Thank you, come again. " )
end ; filter_waves
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191206/b127a3f2/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Screenshot from 2019-12-06 10-08-46.png
Type: image/png
Size: 79674 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191206/b127a3f2/attachment.png>
More information about the ncl-talk
mailing list