[gradsusr] Bias Correct 360days/year files???

Braschi, Andrea andrea.braschi at wur.nl
Wed Jul 25 05:35:02 EDT 2012


Good Morning Grads Users
I am using v2.01.oga.1
I would like to Bias Correct Ensembles (http://ensemblesrt3.dmi.dk/data/A1B/) .tas files
So far I managed (running a Biascorrection script) with some models but not with others 
Interestingly the ones I could not manage are the 360 days/year models (that is my hypothesis for the error i get) 
I noticed that grads cannot overcome this problem (I am dealing with DAILY means)
But i thought to proceed anyway (i am just subtracting the average difference from E-obs file)
The error I get (at file 2 so at model file) while do the bias correction is (in my Bias correction script I assume file 1 is E-obs dataset file 2 is the model.ctl for 20thcentury and file 3 is the model.ctl for 21thcentury)
**********************
NetCDF Error (gancrow, nc_get_vara_double): NetCDF: Index exceeds dimension bound 
Data Request Error: Error for variable 'temp'
	Error occurred at column 1
Syntax Error: Invalid Operand 'tc' not a variable or function name
********************

Temp is defined by me in the .ctl file (here below the ctl file of C4IRCA model that is correctly opened by grads)
-----------------------------------------------------------------------------------
	dset e:\andrea\data\C41\tas\C4IRCA3_A1B_HadCM3Q16_DM_25km_%ch_tas.nc
	chsub 1 3653 1951-1960
	chsub 3654 7305 1961-1970
	chsub 7306 10958 1971-1980
	chsub 10959 14610 1981-1990
	chsub 14611 18623 1991-2000
	chsub 18624 21915 2001-2010
	chsub 21916 25568 2011-2020
	chsub 25569 29220 2021-2030
	chsub 29221 32873 2031-2040
	chsub 32874 36525 2041-2050
	chsub 36526 40178 2051-2060
	chsub 40179 43830 2061-2070
	chsub 43831 47483 2071-2080
	chsub 47484 51135 2081-2090
	chsub 51136 54788 2091-2100
	dtype netcdf
	title C4IRCA3
	undef 1.0E30
	options template
	pdef 174 190 rotll -162.0 39.25 0.22 0.22 -26.16 -20.68
	xdef 464 linear -40.375 0.25
	ydef 201 linear 25.375 0.25
	tdef 54788 linear 0Z01jan1951 1dy
	zdef 1 linear 1 1
	vars 1
	tas=>temp 0 t,z,y,x Mean daily temperature (K)
	endvars
-----------------------------------------------------------------------------------
tc is defined by me into the BiasCorrection scripts here below: ('define Tc=temp-273.15')
-----------------------------------------------------------------------------------
'reinit'
'open E:\andrea\obs\tas\tas_0.25deg_reg.ctl'
'open E:\andrea\data\C41\tas\prova-DMI-C4IRCA3_A1B_DM.ctl'
'open E:\andrea\data\C41\tas\prova-DMI-C4IRCA3_A1B_DM.ctl'
'set vpage 0 11 4 8'
'set grads off'
'set dfile 1'
say 'computing file 1 avg+std.....'

'run E:\andrea\scripts\StatTim tas/100 7671 18628'
'define avgavg1=avgavg'
'define avgstd1=stdavg'
'd avgavg1'
a1=subwrd(result,4)
'd avgstd1'
s1=subwrd(result,4)
'draw title E-OBS observed var: 'a1' +/- 's1
a1

* do second file to compare
'set vpage 0 11 0 4'
'set grads off'
say 'computing file 2 avg+std.....'
'set dfile 2'

* when file 2 is RCM data units for precip are kg.m-2.s-1
* so multiply by 2628000 to get mm/month (avg mo=30.4166 days) or by 86400 to get mm/day
* when file 2 is RCM data units for temperature are K
* so substract 273.15
* do that for whole data set, so set dimensions first

'set t 7306 18263'
'set lon 4.7625'
'set lat 52.375'
*'define pre2=pr*86400'
*'define Tc=tas-273.15'
'define Tc=temp-273.15'

'run E:\andrea\scripts\StatTim Tc 7306 18263'
'define avgavg2=avgavg'
'define avgstd2=stdavg'
'd avgavg2'
a2=subwrd(result,4)
'd avgstd2'
s2=subwrd(result,4)
'draw title RCM-OBS simulated var: 'a2' +/- 's2
'set vpage off'

* give output to screen

say 'avg E-OBS: 'a1' +/- 's1
say 'avg RCM - OBS: 'a2' +/- 's2

* so bias means is difference a1-a2
* so bias stdev is ratio s1/s2

bmean = a1 - a2
smean = a1/a2
bstd  = s1/s2

say 'bias means (diff): 'bmean
say 'bias means (ratio): 'smean
say 'bias std dev (ratio): 'bstd
prompt 'use delta (d) or scaling (s) correction ?  '
pull choice

* do bias correction means
'set dfile 3'
'set lon 4.7625'
'set lat 52.375'
'set t 30682 41639'
* do any transformation when needed (e.g. K to C, etc)
*'define Vc=var1*2628000'
'define Vc=temp-273.15'

if ((choice='d')|(choice='D'))
  say 'doing means-bias correction file 3.....'
* means correction
  'define VcorM= Vc+'bmean
* variance correction
  say 'doing variance-bias correction file 3.....'
  'define VcorMavg=ave(VcorM,t=30682,t=41639)'
  'define delta=VcorM-VcorMavg'
  'define VcorMV=delta*'bstd'+VcorMavg'
else
  say 'doing scaling correction on file 3.....'
* in case of rain, better do a scaling instead of diff correction
  'define VcorMV= Vc*'smean
endif

say 'correction done'
say 'VcorMV is your Variable corrected for bias in Mean and Variance'

'c'
say 'checking...'
say 'original scenario: '
* plot original plus corrected vars to screen
'set vpage 0 11 4 8'
'set grads off'
'set gxout line'
'run E:\andrea\scripts\StatTim Vc 30682 41639'
'd vc'
'draw title RCM-scenA1B simulated original var'
'set vpage 0 11 0 4'
'set grads off'
say 'corrected scenario: '
'run E:\andrea\scripts\StatTim VcorMV 30682 41639'
'd VcorMV'
'draw title RCM-scenA1B simulated bias corrected var'
say 'all done'
--------------------------------------------------------------------------------------------------------------
Here is the Stattim script called into the above one that I have always used for 365,25 days/years models files and always worked so my hypothesis is that cannot be used into 360 day year (3600 time step per decade ) files....
---------------------------------------------------------------------------------------------------------------------------
if (args='')
* if no arguments then ask for variable
  prompt 'which variable ?  '
  pull var
* ask for time interval
  prompt 'begin time ?  '
  pull t1
  prompt 'end time ?  '
  pull t2
else
* get arguments
  var = subwrd(args,1)
  t1  = subwrd(args,2)
  t2  = subwrd(args,3)
endif


'set t 1'
'set lon 4.766667'
'set lat 52.3'
'set z 1'
'define avgavg=ave('%var%',t='%t1%',t='%t2%')'
*'set lon 4.766667'
*'set lat 52.3'
*'define sd=pow('%var%'-avgavg,2)'
*'set x 1'
*'set y 1'
'define stdavg=sqrt(ave(pow('%var%'-avgavg,2),t='%t1%',t='%t2%'))'
'set gxout linefill'
'set lfcols 0 2'
*'set axlim 0.000 0.0001'
*'d avg;avg+std'
*'d avg-std;avg'
*'set gxout line'
*'set ccolor 3'
*'d avg'

*time averaged spatial averages
*'set t 1'
*'define avgavg=ave(avg,t='t1',t='t2')'
'd avgavg'
a1=subwrd(result,4)

*time variance of spatial averages
*'set t '%t1%' '%t2
*'define sd=pow(avg-'a1',2)'
*'set t 1'
*'define stdavg=sqrt(ave(sd,t='t1',t='t2'))'
'd stdavg'
s1=subwrd(result,4)
* time averaged spatial variances
*'set t 1'
*'define avgstd=ave(std,t='t1',t='t2')'
*'d avgstd'
*a2=subwrd(result,4)
*time variance of spatial variances
*'set t '%t1%' '%t2
*'define sd=pow(std-'a2',2)'
*'set t 1'
*'define stdstd=sqrt(ave(sd,t='t1',t='t2'))'
*'d stdstd'
*s2=subwrd(result,4)
*say 'time averaged spatial average+std : 'a1' +/- 's1
*say 'time averaged spatial variance+std: 'a2' +/- 's2	
-----------------------------------------------------------------------------------------------------

Any suggestions?
Thank you

Andrea

Wageningen University 
Msc Student








More information about the gradsusr mailing list