Thanks for reply Dennis.
>What datasets were used?
As in the website (http://www.ncl.ucar.edu/Applications/mjoclivar.shtml), I used http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.html, for example for wind I used daily wind data (http://www.esrl.noaa.gov/psd/cgi-bin/db_search/DBListFiles.pl?did=32&tid=34020&vid=666).
For calculating the anomaly, I just simply modified mjoclivar_2.ncl. What I worried is about accumulating the daily data and then calculated anomaly from accumulated data. If anomaly is calculated from daily data (not from accumulated data) and then daily caclculated anomaly was accumulated, of course, the result will be different if it is calculated from accumulated data.
All plot for OLR is the same as the example, only for wind is different and as we know olr data already accumulated by NOAA, but not the case for wind. Please see the code I used below.
Again, thank you.
Regards,
Joeky
-----the 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 = 19800101 ; start yyyymmdd
ymdLast = 20120303 ; last
yrStrt = ymdStrt/10000
yrLast = ymdLast/10000
lev = 200 ; desired level
nhar = 4 ; number of fourier comp
var = "uwnd" ; name of file
vName = "U" ; 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 ./Reana1/uwnd/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= "./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[0]->lat
lon = f[0]->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.20011001)
ntEnd = ind(yyyymmdd.eq.20020401)
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 <ncl-talk_at_ucar.edu>
Sent: Wednesday, June 6, 2012 11:20 AM
Subject: Re: Wind anomaly, MJO_PC_INDEX and MJO clivar
The real key is the daily anomalies. If they are not correct, none
of the MJO diagnostics are correct. You should look very carefully
at how the anomalies were generated. What datasets were used?
For debugging and testing, use a smaller temporal sample and
plot the data.
On 6/3/12 10:58 PM, juki juki wrote:
> Dear All;
> As my previous email that I would like to plot MJO phase diagram by
> usinghttp://www.ncl.ucar.edu/Applications/Scripts/mjoclivar_15.ncl. To
> see the performance of the code, first I plot same data period as in the
> website. The result is very different from the example
> (http://www.ncl.ucar.edu/Applications/Images/mjoclivar_15_lg.png) in
> which the location of all MJO phases for all months is different from
> that in the website.
> Thus, my MJO_PC_INDEX.nc <http://mjo_pc_index.nc/> may NOT be correct.
> As we know the input of mjoclivar_15.ncl is MJO_PC_INDEX.nc
> <http://MJO_PC_INDEX.nc> which is creted by
> http://www.ncl.ucar.edu/Applications/Scripts/mjoclivar_14.ncl. Moreover,
> the input of
> http://www.ncl.ucar.edu/Applications/Scripts/mjoclivar_14.ncl are:
>
>
> filolr = "olr.day.anomalies.1980-2012.nc
> <http://olr.day.anomalies.1980-2012.nc/>"
> filu200 = "uwnd.day.200.anomalies.1980-2012.nc
> <http://uwnd.day.200.anomalies.1980-2012.nc/>"
> filu850 = "uwnd.day.850.anomalies.1980-2012.nc
>
<http://uwnd.day.850.anomalies.1980-2012.nc/>"
> Among these 3 inputs, it seems the wind anomalies are not correct because
> when I try to calculate Conventional (covariance) univariate EOF
> analysis <http://www.ncl.ucar.edu/Applications/Scripts/mjoclivar_12.ncl>
> fortime span 1995-1999 (as in the website) using
> http://www.ncl.ucar.edu/Applications/Scripts/mjoclivar_12.ncl, the
> result of OLR is close to the website, but for U850 and U200 is
> different and significant different for EOF-3 (enclosed). Thus, I think
> my wind anomalies is not correctly calculated.
> I calculated wind anomalies from daily data as the code below
> (enclosed). I don’t know whether the way to merge the files and then
> calculate anomaly from those files is correct or not. Probably you have
>
suggestion and help for that.
> Thanks you,
> JoeKy
>
>
> _______________________________________________
> 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 Jun 6 03:16:01 2012
This archive was generated by hypermail 2.1.8 : Wed Jun 06 2012 - 15:17:44 MDT
