day2mon regriding
Semyon Grodsky
senya at ATMOS.UMD.EDU
Tue Jan 12 13:19:01 EST 2010
Elisa,
For time averaging or interpolation I use the save_time.gs script (see
attachments).
USAGE:
save_time daily.ctl any_monthly.ctl output_filename
output_filename must be without extension
OUTPUT: ouput_filename.ctl and output_filename.dat
Copy the three *.gs files into your current directory or your
grads_script directory.
Hope it helps.
--Senya
On Fri, 8 Jan 2010, MARIA ELISA SIQUEIRA SILVA wrote:
> Hi,
>
> does anybody know how to construct a gs
> that regrids day or pentad data to monthly data?
>
> thanks
> Elisa
> --
> Maria Elisa Siqueira Silva
> Laboratório de Climatologia e Biogeografia - LCB
> Depto Geografia - FFLCH - USP
> Av. Prof. Lineu Prestes, 338 - Cidade Universitária
> São Paulo - SP
> 05508-900
> (11) 3091-3769
-----------------------------------------------------------
Semyon Grodsky
Computer and Space Science Building (Bldg. #224), Room 2409
Department of Atmospheric and Oceanic Science
University of Maryland
College Park, MD 20742
Phone: 301-405-5330
Fax: 301-314-9482
E-mail: senya at atmos.umd.edu
-------------- next part --------------
function save_time (args)
*'set warn off'
ver = '12/24/2009'
'set warn off'
'last_file'
nf=result
if nf>0
'q file'
res=sublin(result,1)
nfc=subwrd(res,2)
'gxstate1'
dimensions=result
endif
* Get arguments
if (args='')
say '-----------------------------------------------------------------'
say 'This sript changes time decimation '
say 'from INPUT.CTL to that defined by REFERENCE.CTL'
say '-----------------------------------------------------------------'
say 'save_time INPUT.CTL REFERENCE.CTL OUTPUTFILE{no extension} IT a b'
say '-----------------------------------------------------------------'
say 'IT=1 (DEFAULT) --> time averaging'
say 'use IT=1 if the new time step exceeds the original time step'
say ''
say 'IT=0 --> time interpolation'
say 'use IT=0 if the new time step is less than the original time step'
say '-----------------------------------------------------------------'
say 'output=a*var+b, DEFAULT: a=1, b=0'
say '-----------------------------------------------------------------'
say 'ver 'ver' by senya at atmos.umd.edu'
say '-----------------------------------------------------------------'
return
else
inputfile = subwrd(args,1)
referencefile = subwrd(args,2)
outputfile=subwrd(args,3)
itime=subwrd(args,4)
a=subwrd(args,5)
b=subwrd(args,6)
endif
***********************************************************
* referencefile='/mnt/data/senya/pacific/sst/sst_month.ctl'
* referencefile='/mnt/data/senya/pacific/sst/year.ctl'
***********************************************************
if a=''
a=1
endif
if b=''
b=0
endif
if itime=''
itime=1
endif
'get_path 'outputfile
path=result
'get_filename 'outputfile
fnameout=result
say '********************************************************'
say 'SAVE_TIME ver: 'ver
say '--------------------------------------------------------'
say 'INPUT CTL : 'inputfile
say 'REFERENCE CTL: 'referencefile
say 'OUTPUT CTL : 'outputfile%'.ctl'
say '********************************************************'
**************************************if 1=2
* 'reinit'
'open 'inputfile
'set dfile 'nf+1
'get_undef'
undef=result
'open 'referencefile
res=sublin(result,2)
if subwrd(res,1)='Open'
say result
say 'Check referencefile'
return
endif
* 'q files'
* say result
* pull dummy
'q file'
res=sublin(result,5)
xl=subwrd(res,3)
yl=subwrd(res,6)
zl=subwrd(res,9)
res=sublin(result,6)
nv=subwrd(res,5)
'set x 1 'xl
'set y 1 'yl
'set z 1 'zl
'set t 1 last'
'q time'
res=subwrd(result,3)
time1=substr(res,6,7)
res=subwrd(result,5)
time2=substr(res,6,7)
'set dfile 'nf+2
'q dims'
res=sublin(result,5)
tm1=subwrd(res,11)
tm2=subwrd(res,13)
tm1=math_int(tm1)
tm2=math_int(tm2)
say tm1' 'tm2
*pull dummy
'set gxout fwrite'
'set fwrite 'outputfile'.dat'
i=tm1
while(i<=tm2)
'set dfile 'nf+2
'set t 'i
'q time'
t1=subwrd(result,3)
if i=tm1
timestart=t1
endif
if i=tm2
timeend=t1
endif
i1=i+1
'set t 'i1
'q time'
t2=subwrd(result,3)
'set t 'i
'set dfile 'nf+1
'q dims'
res=sublin(result,5)
qt=subwrd(res,9)
iqt1=math_int(qt)
wt1=1-(qt-iqt1)
iqt2=iqt1+1
wt2=1-wt1
****************
'set t 1'
****************
iv=1
while(iv<=nv)
'q file'
res=sublin(result,6+iv)
var=subwrd(res,1)
zl=subwrd(res,2)
if zl=0
zl=1
endif
'set z 1 'zl
if i=tm1
say var' levels='zl
endif
*case time averaging
if itime=1
say t1 ' ' t2
'varr=ave('a'*'var',time='t1',time='t2')+'b
endif
if itime=0
say t1 '<== 'iqt1' 'iqt2' 'wt1' 'wt2
'varr='var'(t='iqt1')*'wt1'+'var'(t='iqt2')*'wt2
'varr='a'*varr+'b
endif
'varr=const(varr, 'undef', -u)'
*'d varr'
*----------------------------
iz=1
while (iz<=zl)
'set z 'iz
'd varr'
iz=iz+1
endwhile
'set z 1 'zl
*----------------------------
iv=iv+1
endwhile
i=i+1
endwhile
'disable fwrite'
'set gxout contour'
say '----------------------------'
say timestart' 'timeend
nm=tm2-tm1+1
say 'step1=' tm1 ' step2=' tm2 ' timesteps='nm
*find new timestep
i=1
while(i<1000)
res=read(referencefile)
outlin=sublin(res,2)
par=subwrd(outlin,1)
if (par='TDEF'|par='tdef'|par='Tdef')
newtimestep=subwrd(outlin,5)
i=1000
endif
i=i+1
endwhile
**********************************
*save CTL
*************endif
*****************
*nm=0
*********************
i1=0
i2=0
i=1
status=0
while(status = 0)
res=read(inputfile)
status=sublin(res,1)
outlin=sublin(res,2)
par=subwrd(outlin,1)
if (par='OPTIONS'|par='options'|par='Options')
say '==> OPTIONS LINE is REMOVED'
outlin=' '
endif
if (par='FORMAT'|par='format'|par='Format')
say '==> FORMAT LINE is REMOVED'
outlin=' '
endif
if (par='DTYPE'|par='dtype'|par='Dtype')
say '==> DTYPE LINE is REMOVED'
outlin=' '
endif
if (par='UNPACK'|par='unpack'|par='Unpack')
say '==> UNPACK LINE is REMOVED'
outlin=' '
endif
if (par='DSET'|par='dset'|par='Dset')
outlin='DSET ^'fnameout'.dat'
endif
if (par='TDEF'|par='tdef'|par='Tdef')
outlin='TDEF ' math_int(nm) ' LINEAR ' timestart ' 'newtimestep
endif
if (par='ENDVARS'|par='endvars'|par='Endvars')
outlin='ENDVARS'
endif
if (par='VARS'|par='vars'|par='Vars')
i1=i+1
ni=subwrd(outlin,2)
i2=i+ni
endif
if i=1
dum=write(outputfile'.ctl',outlin)
else
if status=0 & (i-i1)*(i-i2)>0
dum=write(outputfile'.ctl',outlin, append)
endif
if status=0 & (i-i1)*(i-i2)<=0
******change vars format to float **********
outlin_new=''
j=1
while subwrd(outlin,j)!=''
if j!=3
outlin_new = outlin_new%subwrd(outlin,j)%' '
else
outlin_new = outlin_new%'99 '
endif
j=j+1
endwhile
dum=write(outputfile'.ctl',outlin_new, append)
endif
***********************************************
endif
i=i+1
endwhile
say '********************************************************'
say 'SAVE_TIME ver: 'ver
say '--------------------------------------------------------'
say 'INPUT CTL : 'inputfile
say 'REFERENCE CTL: 'referencefile
say 'OUTPUT CTL : 'outputfile%'.ctl'
say '********************************************************'
'close 'nf+2
'close 'nf+1
if nf>0
'set dfile 'nfc
'gxstate1 'dimensions
endif
-------------- next part --------------
******************************************************
function gxstate(args)
*
* Get or restore dimensions
* get dimension USAGE: gxstate1
* restore dimension USAGE: gxstate1 dimensions
* dimeansions is the output of get dimension by gxstate1()
******************
if args=''
'query dim'
dinf = result
lx = sublin(dinf,2)
ly = sublin(dinf,3)
lz = sublin(dinf,4)
lt = sublin(dinf,5)
if ( subwrd(lx,7) = 'to')
lons = subwrd(lx,6)
lone = subwrd(lx,8)
xs = subwrd(lx,11)
xe = subwrd(lx,13)
else
lons = subwrd(lx,6)
lone = subwrd(lx,6)
xs = subwrd(lx,9)
xe = subwrd(lx,9)
endif
if ( subwrd(ly,7) = 'to')
lats = subwrd(ly,6)
late = subwrd(ly,8)
ys = subwrd(ly,11)
ye = subwrd(ly,13)
else
lats = subwrd(ly,6)
late = subwrd(ly,6)
ys = subwrd(ly,9)
ye = subwrd(ly,9)
endif
if ( subwrd(lz,7) = 'to')
levs = subwrd(lz,6)
leve = subwrd(lz,8)
zs = subwrd(lz,11)
ze = subwrd(lz,13)
else
levs = subwrd(lz,6)
leve = subwrd(lz,6)
zs = subwrd(lz,9)
ze = subwrd(lz,9)
endif
if ( subwrd(lt,7) = 'to')
tims = subwrd(lt,6)
time = subwrd(lt,8)
ts = subwrd(lt,11)
te = subwrd(lt,13)
else
tims = subwrd(lt,6)
time = subwrd(lt,6)
ts = subwrd(lt,9)
te = subwrd(lt,9)
endif
return(lons' 'lone' 'xs' 'xe' 'lats' 'late' 'ys' 'ye' 'levs' 'leve' 'zs' 'ze' 'tims' 'time' 'ts' 'te)
********end get dimension***********
else
'set lon 'subwrd(args,1)' 'subwrd(args,2)
'set lat 'subwrd(args,5)' 'subwrd(args,6)
'set lev 'subwrd(args,9)' 'subwrd(args,10)
*'set time 'subwrd(args,13)' 'subwrd(args,14)
'set t 'subwrd(args,15)' 'subwrd(args,16)
endif
*****************************************************
-------------- next part --------------
function last_file ()
'q files'
f=result
if subwrd(result,1)='No'
return(0)
else
ln='qqq'
i=1
while (ln != '')
ln=sublin(f,i)
if subwrd(ln,1)='File'
nf=subwrd(ln,2)
endif
i=i+1
endwhile
return(nf)
endif
More information about the gradsusr
mailing list