[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