Re: MJO clivar

From: juki juki <juky_emc2_at_nyahnyahspammersnyahnyah>
Date: Wed May 30 2012 - 19:31:00 MDT

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