Thanks for reply Dennis;
I tried to modify Example 2 on the MJO Application page. The code is given below. I hope the way to merge many files being correct. I have some an error, particularly when writing the output to a file. The errors are:
fatal:Dimension sizes on right hand side of assignment do not match
dimension sizes of left hand side
fatal:Execute: Error occurred at or near line 154 in mjoclivar_2uwnd.ncl (the following line: fnc->lat = (/lat/))
Probably you know the reason..
Thank you
JuKY
---------------------------code
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
begin
;-------------------------------------------------------------
; User specifications
;-------------------------------------------------------------
NC = True ; create netCDF?
PLOT = False ; sample plots?
ymdStrt = 20010101 ; start yyyymmdd
ymdLast = 20111231 ; last
yrStrt = ymdStrt/10000
yrLast = ymdLast/10000
lev = 850 ; desired level
nhar = 4 ; number of fourier comp
var = "uwnd" ; name of file
vName = "UWND" ; name for plots
fili = var+".day.mean.nc" ; input file
;========================
; get list of all files and open as "one big file"
;========================
all_files = systemfunc ("ls e:/Rdata/MJO_cal/Data/uwnd*.nc")
f = addfiles (all_files, "r") ; note the "s" of addfile
print(all_files)
;========================
; choose how files are combined and read in variable across files
;========================
ListSetType (f, "cat") ; concatenate or "merge" (default)
if (NC) then
diro= "e:/RData/MJO_cal/Anom/" ; output dir
filo= var+".day."+lev+".anomalies."+yrStrt+"-"+yrLast+".nc" ; output file
end if
if (PLOT) then
wksType = "ps"
wksName = "mjoclivar" ; "mjo."+yrStrt+"_"+yrLast
end if
;***********************************************************
; Read user specified time and create required yyyyddd
;***********************************************************
; f = addfile (diri+fili , "r")
time = f[:]->time ; all times on file
ymd = cd_calendar(time, -2) ; yyyymmdd
iStrt = ind(ymd.eq.ymdStrt) ; index start
iLast = ind(ymd.eq.ymdLast) ; index last
delete(time)
delete(ymd)
;***********************************************************
; Read user specified time and create required yyyyddd
;***********************************************************
time = f[:]->time(iStrt:iLast) ; time:units = "hours since"
TIME = cd_calendar(time, 0) ; type float
year = floattointeger( TIME(:,0) )
month = floattointeger( TIME(:,1) )
day = floattointeger( TIME(:,2) )
ddd = day_of_year(year, month, day)
yyyyddd = year*1000 + ddd ; needed for input
;***********************************************************
; Read data: short2flt
;***********************************************************
x = short2flt(f[:]->$var$(iStrt:iLast,{lev},:,:)) ; convert to float
printVarSummary(x)
;***********************************************************
; Compute daily climatology: raw and then 'smoothed'
;***********************************************************
xClmDay = clmDayTLL(x, yyyyddd) ; daily climatology at each grid point
printVarSummary(xClmDay)
;***********************************************************
; Compute smoothed daily climatology using 'nhar' harmonics
;***********************************************************
xClmDay_sm = smthClmDayTLL(xClmDay, nhar)
printVarSummary(xClmDay_sm)
;***********************************************************
; Compute daily anomalies using raw and smoothed climatologies
;***********************************************************
xAnom = calcDayAnomTLL (x, yyyyddd, xClmDay)
printVarSummary(xAnom)
printMinMax(xAnom, True)
xAnom_sm = calcDayAnomTLL (x, yyyyddd, xClmDay_sm)
xAnom_sm_at_long_name = "Anomalies from Smooth Daily Climatology"
printVarSummary(xAnom_sm)
printMinMax(xAnom_sm, True)
delete( x ) ; no longer needed
;***********************************************************
; Create netCDF: convenience use 'simple' method
;***********************************************************
dimx = dimsizes(xAnom)
ntim = dimx(0)
nlat = dimx(1)
mlon = dimx(2)
if (NC) then
lat = f[:]->lat
lon = f[:]->lon
system("/bin/rm -f "+diro+filo) ; rm any pre-exist file, if any
fnc = addfile (diro+filo, "c")
filAtt = 0
filAtt_at_title = vName+": Daily Anomalies: "+yrStrt+"-"+yrLast
filAtt_at_source_file = fili
filAtt_at_creation_date = systemfunc("date")
fileattdef( fnc, filAtt ) ; copy file attributes
setfileoption(fnc,"DefineMode",True)
varNC = vName+"_anom"
varNC_sm = vName+"_anom_sm"
dimNames = (/"time", "lat", "lon"/)
dimSizes = (/ -1 , nlat, mlon/)
dimUnlim = (/ True , False, False/)
filedimdef(fnc,dimNames,dimSizes,dimUnlim)
filevardef(fnc, "time" ,typeof(time),getvardims(time))
filevardef(fnc, "lat" ,typeof(lat) ,getvardims(lat))
filevardef(fnc, "lon" ,typeof(lon) ,getvardims(lon))
filevardef(fnc, varNC ,typeof(xAnom) ,getvardims(xAnom))
filevardef(fnc, varNC_sm,typeof(xAnom) ,getvardims(xAnom))
filevarattdef(fnc,"time" ,time) ; copy time attributes
filevarattdef(fnc,"lat" ,lat) ; copy lat attributes
filevarattdef(fnc,"lon" ,lon) ; copy lon attributes
filevarattdef(fnc,varNC ,xAnom)
filevarattdef(fnc,varNC_sm,xAnom_sm)
fnc->time = (/time/)
fnc->lat = (/lat/)
fnc->lon = (/lon/)
fnc->$varNC$ = (/xAnom /)
fnc->$varNC_sm$ = (/xAnom_sm/)
end if
;************************************************
; plotting parameters
;************************************************
if (PLOT) then
LAT = (/ 60, 45, 5, -5, -45, 60 /)
LON = (/270, 30, 90, 90, 180, 0 /)
nPts = dimsizes( LAT )
plot = new ( nPts, graphic)
data = new ( (/2,366/), typeof(xClmDay), getFillValue(xClmDay))
wks = gsn_open_wks (wksType,wksName)
res = True ; plot mods desired
res_at_gsnDraw = False
res_at_gsnFrame = False
res_at_trXMinF = 1
res_at_trXMaxF = 366
;res_at_tiMainString = ""
res_at_xyLineThicknesses = (/1.0, 2.0/) ; make 2nd lines thicker
res_at_xyLineColors = (/"blue","red"/) ; change line color
res_at_xyMonoDashPattern = True ; all solid
do np=0,nPts-1
data(0,:) = xClmDay(:,{LAT(np)},{LON(np)})
data(1,:) = xClmDay_sm(:,{LAT(np)},{LON(np)})
res_at_gsnCenterString = "lat="+LAT(np)+" lon="+LON(np)
plot(np) = gsn_csm_y (wks,data,res) ; create plot
end do
resP = True ; modify the panel plot
resP_at_txString = vName+": Raw/Smooth Daily Climatology: "+yrStrt+"-"+yrLast
resP_at_gsnMaximize = True
resP_at_gsnPaperOrientation = "portrait"
gsn_panel(wks,plot,(/(nPts/2),2/),resP) ; now draw as one plot
;==========
; Plot anomalies for an arbitrarily selected near equatorial location
; Time: Oct 1, 1996 to April 1,1997 [arbitrary selection]
;==========
LATX = 0
LONX = 90
yyyymmdd = cd_calendar(time, -2)
;;yrfrac = yyyymmdd_to_yyyyfrac (yyyymmdd, 0)
;;delete(yrfrac_at_long_name)
xAnom_at_long_name = "Anomalies from Raw" ; short labels for plot
xAnom_sm_at_long_name = "Anomalies from Smooth"
ntBegin = ind(yyyymmdd.eq.19961001)
ntEnd = ind(yyyymmdd.eq.19970401)
monthLabels = (/1,4,7,10/)
monNam = (/"Jan","Feb","Mar","Apr","May","Jun" \
,"Jul","Aug","Sep","Oct","Nov","Dec" /)
xVal = new(ntim, typeof(xAnom&time) , "No_FillValue") ; bigger than
xLab = new(ntim, "string", "No_FillValue") ; needed
xValm = new(ntim, typeof(xAnom&time) , "No_FillValue") ; bigger than
ntm = -1
cr = inttochar(10) ; carriage return
do nt=ntBegin,ntEnd
if (day(nt).eq.1) then
ntm = ntm + 1
xVal(ntm) = xAnom&time(nt)
xLab(ntm) = monNam(month(nt)-1)
if (month(nt).eq.1) then
xLab(ntm) = xLab(ntm) + cr +sprinti("%0.4i", year(nt))
end if
end if
end do
rxy = True
rxy_at_gsnDraw = False
rxy_at_gsnFrame = False
rxy_at_gsnYRefLine = 0.0 ; create a reference line
rxy_at_gsnAboveYRefLineColor = "red" ; above ref line fill red
rxy_at_gsnBelowYRefLineColor = "blue" ; below ref line fill blue
rxy_at_xyLineThicknessF = 2.0
rxy_at_vpHeightF = 0.4 ; resize
rxy_at_vpWidthF = 0.8
rxy_at_tmXBMode = "Explicit"
rxy_at_tmXBValues = xVal(0:ntm)
rxy_at_tmXBLabels = xLab(0:ntm)
plot(0) = gsn_csm_xy (wks,time(ntBegin:ntEnd) \
,xAnom(ntBegin:ntEnd,{0},{90}),rxy)
plot(1) = gsn_csm_xy (wks,time(ntBegin:ntEnd) \
,xAnom_sm(ntBegin:ntEnd,{0},{90}),rxy)
resP_at_txString = vName+": Anomalies: (0,90E)"
gsn_panel(wks,plot(0:1),(/2,1/),resP)
end if
end
________________________________
From: Dennis Shea <shea_at_ucar.edu>
To: juki juki <juky_emc2_at_yahoo.com>
Cc: "ncl-talk_at_ucar.edu" <ncl-talk_at_ucar.edu>
Sent: Thursday, May 17, 2012 10:16 AM
Subject: Re: [ncl-talk] MJO clivar
[1] There is no function to calculate the propagation speed
of the MJO.
[2] "...compare the result with another method/result".
Not sure what other method you are referring to.
It would be a surprise if the other method differed in any
substantive way with that in MJO_CLIVAR.
[3] There are assorted examples on the Application pages
that can be used as a basis for what you need. In paarticular,
Example 2 on the MJO Application page
http://www.ncl.ucar.edu/Applications/mjoclivar.shtml
Please read the documentation of the assorted functions
and syntax ( {...} )used.
Many climate/anomaly functions are at:
http://www.ncl.ucar.edu/Document/Functions/climo.shtml
[4] The following is an excerpt and modification of example 2,
Note the use of addfiles and the {...}/ $...$ syntax
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
var = "uwnd"
ymdStrt = 19800101 ; start yyyymmdd
ymdLast = 20111231 ; last
yrStrt = ymdStrt/10000
yrLast = ymdLast/10000
diri = "./"
fils = systemfunc("ls " + diri + "uwnd*")
print(fils) ; full path to files
latS = -10 ; if desired, specify a region and level
latN = 10
lonL = 90
lonR = 270
lev = 200 ; desired level
f = addfiles( diri+fils, "r")
time = f[:]->time ; all times on files
ymd = cd_calendar(time, -2) ; yyyymmdd
iStrt = ind(ymd.eq.ymdStrt) ; index start
iLast = ind(ymd.eq.ymdLast) ; index last
delete([/time, ymd/])
;***********************************************************
; Read user specified time; create required yyyyddd
;***********************************************************
time = f->time(iStrt:iLast) ; time:units = "hours since"
TIME = cd_calendar(time, 0) ; type float
year = floattointeger( TIME(:,0) )
month = floattointeger( TIME(:,1) )
day = floattointeger( TIME(:,2) )
ddd = day_of_year(year, month, day)
yyyyddd = year*1000 + ddd ; needed for input
x = short2flt(f[:]->$var$(iStrt:iLast,{lev},{latS:latN},{lonL:lonR}))
printVarSummary(x) ; (time,lat,lon)
xClm = clmDayTLL(x) ; (366,lat,lon)
xAnom= calcDayAnomTLL ( x, xClm)
printVarSummary(xClm)
prntVarSummary(x, uClm)
The 'xAnom' would be uwnd.day.200.anomalies.1980-2011.nc
or
var+".day."+lev+".anomalies."+yrStrt+"-"+yrLast+".nc"
The netCDF can be created analogous to Example 2
On 5/16/12 5:44 AM, juki juki wrote:
> Dear All;
>
> I need to make the MJO clivar from 1980-2010 data (OLR). I read the
> following link to do it
> (http://www.ncl.ucar.edu/Applications/Scripts/mjoclivar_14.ncl).
> However, there are 3 data type in it i.e.,
>
> filolr = "olr.day.anomalies.1980-2005.nc <http://olr.day.anomalies.1980-2005.nc/>"
> filu200 = "uwnd.day.200.anomalies.1980-2005.nc <http://uwnd.day.200.anomalies.1980-2005.nc/>"
> filu850 = "uwnd.day.850.anomalies.1980-2005.nc <http://uwnd.day.850.anomalies.1980-2005.nc/>"
>
>
> For filolr I can produce it because the daily olr data already in one file, but for"uwnd.day.200.anomalies.1980-2005.nc", I think uwnd daily data is one file/year (http://www.esrl.noaa.gov/psd/cgi-bin/db_search/DBListFiles.pl?did=59&tid=32778&vid=1524). Probably any body has suggestion to read many file sequently and then produce"uwnd.day.200.anomalies.1980-2005.nc".
>
> Another question, is it any ncl tool to calculate the propagation speed of MJO ? Currently
> I am working with the tracking of MJO propagation. It is worthwhile to compare the result with another method/result.
>
>
> Thanks for sharing,
>
>
> JuKY
>
>
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed May 30 19:31:10 2012
This archive was generated by hypermail 2.1.8 : Wed Jun 06 2012 - 15:17:44 MDT