average.
Charles Seman
Charles.Seman at NOAA.GOV
Mon Oct 3 18:58:26 EDT 2005
Dear Diego Medina,
Please find below a copy of a message (copied from a "forward" email
message window) which I sent to the GrADS Listserv for a similar problem.
Perhaps you could try something like (note, I haven't tested the
following code):
month_days = '31 28 31 30 31 30 31 31 30 31 30 31'
*
* loop over months...
*
m=1
while ( m <= 12 )
*
* define number of days for this month
*
ndays = subwrd(month_days,m)
*
* open first day of month file and use it to initialize a time sum
variable...
*
ctl_file = '1990'm'01.ctl' * assuming you are using a GrADS ctl file
with a name like this; correct if not so...
'open 'ctl_file
'query file'
say result
say
rec5 = sublin(result,5) ; rec6 = sublin(result,6)
xsiz = subwrd(rec5,3) ; ysiz = subwrd(rec5,6) ; zsiz = subwrd(rec5,9)
; tsiz = subwrd(rec5,12)
nvar = subwrd(rec6,5)
say ' nvar 'nvar
say ' xsiz 'xsiz
say ' ysiz 'ysiz
say ' zsiz 'zsiz
say ' tsiz 'tsiz
*
'set x 1 ' xsiz
'set y 1 ' ysiz
'set z 1 ' zsiz
'set t 1' * wildcard for time sum and average...
'define usum = ave(u,t=1,t='tsiz')' * if you have more than one time
level per daily file...
'define vsum = ave(v,t=1,t='tsiz')' * if you have more than one time
level per daily file...
'define usum = u' * if you have one time level per daily file...
'define vsum = v' * if you have one time level per daily file...
*
* close file for this day...
*
'close 1'
*
* update "sums" for each of the rest of the days for this month...
*
n=2
while ( n <= ndays )
*
* define ctl file in terms of day index "n"
* (assuming you are using a GrADS ctl file with a name like this;
correct if not so...)
*
if( n < 10 )
ctl_file = '1990'm'0'n'.ctl'
else
ctl_file = '1990'm''n'.ctl'
endif
'open 'ctl_file
'query file'
say result
say
rec5 = sublin(result,5) ; rec6 = sublin(result,6)
xsiz = subwrd(rec5,3) ; ysiz = subwrd(rec5,6) ; zsiz =
subwrd(rec5,9) ; tsiz = subwrd(rec5,12)
nvar = subwrd(rec6,5)
say ' nvar 'nvar
say ' xsiz 'xsiz
say ' ysiz 'ysiz
say ' zsiz 'zsiz
say ' tsiz 'tsiz
*
'set x 1 ' xsiz
'set y 1 ' ysiz
'set z 1 ' zsiz
'set t 1' * wildcard for time sum and average...
'usum = usum + ave(u,t=1,t='tsiz')' * if you have more than one
time level per daily file...
'vsum = vsum + ave(v,t=1,t='tsiz')' * if you have more than one
time level per daily file...
'usum = usum + u' * if you have one time level per daily file...
'vsum = vsum + v' * if you have one time level per daily file...
*
* close file for this day...
*
'close 1'
n=n+1
endwhile
*
* calculate the average of the daily sums for this month...
*
'define uavg = usum/'ndays
'define vavg = vsum/'ndays
*
* display data: (x,y) plots over the vertical levels shown below...
*
zz=1
while ( zz <= zsiz )
'set x 1 ...'
'set y 1 ...'
'set z 'zz
'set t 1' * wildcard for time average...
'set gxout shaded'
'd uavg'
'd vavg'
zz=zz+1
endwhile
*
* finished for this month...
*
'undefine usum' ; 'undefine uavg'
'undefine vsum' ; 'undefine vavg'
*
* go to next month...
*
m=m+1
endwhile
I hope this helps; please let me know if you encounter any problems
with (or have any questions about) the above code, OK?
Thanks,
Chuck
diego wrote:
>Dear grads User,
>
>I have a mean diary for u and v for sudamerica in
>files :
>1990101.grb
>1990102.grb
>1990103.grb
>1990104.grb
>1990105.grb
>...
>...
>1990131.grb
>
>1990201.grb
>1990202.grb
>...
>...
>
>I want to obtain monthly averages of u and v.
>can you help me.
>
>Thanks
>
>
>
>
>
>______________________________________________
>Renovamos el Correo Yahoo!
>Nuevos servicios, más seguridad
>http://correo.yahoo.es
>
>
------------------------------------------------------------------------------------
-------- Original Message --------
Subject: Re: average in reanalysis data
Date: Fri, 18 Mar 2005 18:22:46 -0500
From: Charles Seman <Charles.Seman at noaa.gov>
To: GRADSUSR at LIST.CINECA.IT
References: <bd4c9dd96376.42385be5 at unal.edu.co>
Luis,
Perhaps you could try calculating an average of 30 years of averaged
precipitation data as follows:
ts = 1961
tf = 1990
nyear = tf-ts+1
*
* initialize "pravgsum" for the first year...
*
tt=ts
*
* define reanalysis file in terms of year index "tt"
*
reanalysis_file = ... * must be defined
' open 'reanalysis_file
'query file'
say result
say ' '
rec5 = sublin(result,5) ; rec6 = sublin(result,6)
xsiz = subwrd(rec5,3) ; ysiz = subwrd(rec5,6) ; zsiz = subwrd(rec5,9) ;
tsiz = subwrd(rec5,12)
nvar = subwrd(rec6,5)
say ' nvar 'nvar
say ' xsiz 'xsiz
say ' ysiz 'ysiz
say ' zsiz 'zsiz
say ' tsiz 'tsiz
*
'set lon ...' * must be defined
'set lat ...' * must be defined
'set lev ...' * must be defined
'set t 1'
'define pravg = ave('precip_variable',t=1,t='tsiz')' * must be defined
*
* start sum of yearly precip averages...
*
'pravgsum = pravg'
'undefine pravg'
*
* close reanalysis file for this year...
*
'close 1'
*
* update "pravgsum" for each of the rest of the years...
*
tt=tt+1
while ( tt <= tf )
*
* define reanalysis file in terms of year index "tt"
*
reanalysis_file = ... * must be defined
'open 'reanalysis_file
'set lon ...' * must be defined
'set lat ...' * must be defined
'set lev ...' * must be defined
'set t 1'
'define pravg = ave('precip_variable',t=1,t='tsiz')' * must be defined
*
* update "pravgsum" using "pravg" for this year...
*
'pravgsum = pravgsum + pravg'
'undefine pravg'
*
* close reanalysis file for this year...
*
'close 1'
tt=tt+1
endwhile
*
* calculate the average of the "nyear" yearly averages...
*
'define pravg = pravgsum/'nyear
'd pravg'
....
I hope this helps... Please let me know if you find any problems with
the above technique, OK?
Thanks,
Chuck
Luis Fernando Montana R. wrote:
> Hello all,
>
> I'm using data of the reanalysis for period 1961-1990. I have 30
> files, one per year and need to calculate the average of precipitation
> for that period.
>
> There is some simplified form to do it?
>
> Thanks in advance,
>
> Luis Fernando Montana Roa
> Posgrado de Meteorologia
> Universidad Nacional de Colombia
> Bogota - Colombia - South America
>
--
********************************************************************
Charles Seman Charles.Seman at noaa.gov
U.S. Department of Commerce / NOAA / OAR
Geophysical Fluid Dynamics Laboratory voice: (609) 452-6547
201 Forrestal Road fax: (609) 987-5063
Princeton, NJ 08542 http://www.gfdl.noaa.gov/~cjs/
********************************************************************
"The contents of this message are mine personally and do not reflect
any position of the Government or NOAA."
------------------------------------------------------------------------------------
--
********************************************************************
Charles Seman charles.seman at noaa.gov
U.S. Department of Commerce / NOAA / OAR
Geophysical Fluid Dynamics Laboratory voice: (609) 452-6547
201 Forrestal Road fax: (609) 987-5063
Princeton, NJ 08540-6649 http://www.gfdl.noaa.gov/~cjs/
********************************************************************
"The contents of this message are mine personally and do not reflect
any position of the Government or NOAA."
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://gradsusr.org/pipermail/gradsusr/attachments/20051003/88c2e024/attachment.html
More information about the gradsusr
mailing list