********************************************************************************
*                                                                              *
*  The GrADS script performs calculation of the following verifications scores *
*  given analysis and prediction results available in GrADS form:              *
* RMSE =                                                                       *
* TEND CORR =                                                                  *
* MEAN ABS ERR =                                                               *
* BIAS =                                                                       *
* S1 =                                                                         *
*  Each time the calculations are performed for a single                       *
*  surface. Title , date, "lev" ,"lat1,lat2,lon1,lon2"                         *
*  parametrs  have to be determined  manually.                                 *
*  Results are directed to outver.txt file in ascii form                       *
********************************************************************************
*  History: Simon Krichak, Melina Dayan, Tel Aviv University                   *
*  (shimon@cyclone.tau.ac.il), June 2004.                                      *
********************************************************************************

h=24
'set lev 1000'

dt=5
dir=120
dtt=dt+1
dttt=dt+2

while(h<72)
* open GrADS file with the forecast result for "h" hours
'open /'2004%dir%dt'/domain1:'2004%dir%dt%_00'.ctl'
* open GrADS file with the initial data for the fcst
'open /'2004%dir%dt'/regrid_domain1:'2004%dir%dt%_00'.ctl'
* open GrADS file with the first future date
'open /'2004%dir%dtt'/regrid_domain1:'2004%dir%dtt%_00'.ctl'
* open GrADS file with the second future date
'open /'2004%dir%dttt'/regrid_domain1:'2004%dir%dttt%_00'.ctl'

'set lev 1000'

'define fcst = pslv.1(t='h/3+1')'
'define init = pslv.2(t=1)/100'

if (h<48)
'define futr = pslv.3(t=1)/100'
else
'define futr = pslv.4(t=1)/100'
endif

lat1 = 30
lat2 = 35

lon1 = 33
lon2 = 38
lat2m1 = lat2 - 1
lon2m1 = lon2 - 1

'set lat 'lat1' 'lat2'' 
'set lon 'lon1' 'lon2''

#' Define real  time variation '
'define c1 = futr - init'

#' Define predicted  time variation '
'define c2 = fcst - init'

#' Define forecast error'
'define c3 = fcst - futr'

# (futr - init)**2
'define s6 =  c1*c1'

# (fcst - init)**2
'define s7 =  c2*c2'

# (fcst - futr)**2
'define s8 =  c3*c3'

# (fcst - init) * (futr - init)
'define s9 =  c1*c2'
 
# mean real tendency
' a =  aave(c1,lon='lon1',lon='lon2',lat='lat1',lat='lat2')'

# mean fcsted tendency
' b =  aave(c2,lon='lon1',lon='lon2',lat='lat1',lat='lat2')'

# mean abs (real  time variation i.e. tendency)
' sma =  aave(abs(c1),lon='lon1',lon='lon2',lat='lat1',lat='lat2')'

# mean abs (predicted  time variation i.e. tendency)
' smb =  aave(abs(c2),lon='lon1',lon='lon2',lat='lat1',lat='lat2')'

# mean abs (forecast error)
' smer = aave(abs(c3),lon='lon1',lon='lon2',lat='lat1',lat='lat2')'

' bias = aave(c3,lon='lon1',lon='lon2',lat='lat1',lat='lat2')'

'  ss6 = aave(s6,lon='lon1',lon='lon2',lat='lat1',lat='lat2')'
'  ss7 = aave(s7,lon='lon1',lon='lon2',lat='lat1',lat='lat2')'
'  ss8 = aave(s8,lon='lon1',lon='lon2',lat='lat1',lat='lat2')'
'  ss9 = aave(s9,lon='lon1',lon='lon2',lat='lat1',lat='lat2')'
'  ska = sqrt(ss6)'
'  skb = sqrt(ss7)'
'  sker = sqrt(ss8)'
'  esk = sker/ska'
  
#____________________
# Relative abs error
' epr =  smer/sma'
#____________________

' rsk = ss9/(sqrt(ss6*ss7))'
' al = skb/ska'

# _________________________________
' define cc1 = init - fcst - a'
' define cc2 = futr - fcst - b' 
' define ss1 = cc1*cc2'
' define ss2 = cc1*cc1'
' define ss3 = cc2*cc2'
' s1 = aave(ss1,lon='lon1',lon='lon2',lat='lat1',lat='lat2')'
' s2 = aave(ss2,lon='lon1',lon='lon2',lat='lat1',lat='lat2')'
' s3 = aave(ss3,lon='lon1',lon='lon2',lat='lat1',lat='lat2')'
' rpr = s1/sqrt(s2*s3)'
#' if (s2 = 0)' 
#' s2 = 0.001'
#' endif'
#' if (s3 = 0)'
#' s3 = 0.001'
#' endif'


'set lat 'lat1' 'lat2m1'' 
'set lon 'lon1' 'lon2m1''

* predicted gradients
' define a = cdiff(fcst,x)' 
' define b = cdiff(fcst,y)'
*real future gradients
' define c = cdiff(futr,x)'
' define d = cdiff(futr,y)'

* abs diff between predicted and future real gradients
' define a1 = abs(a-c)'
' define b1 = abs(b-d)'

* abs -x gradients  
' define c1 = abs(a)'
' define c11 = abs(c)'

* abs -y gradients  
' define c2 = abs(b)'
' define c21 = abs(d)'

* c111 negative if c11 > c1
' define c111 = c1 - c11'
* c112 negative if c11 < c1
' define c112 = c11 - c1'

* c122 = UNDEF if c111 < 0
' define c122 =  maskout(c1,c1 - c11) '
* newc122 = c122 but zero (0) instead of UNDEF values
' define newc122 = const(c1,0,-u)'

* c123 = UNDEF if c112 < 0
' define c123 =  maskout(c11,c11 - c1)'
* newc123 = c123 but zero (0) instead of UNDEF values
' define newc123 = const(c11,0,-u)'

* c211 negative if c21 > c2
' define c211 = c2 - c21'
* c212 negative if c2 > c21
' define c212 = c21 - c2'

* c222 = UNDEF if c211 < 0
' define c222 = maskout(c2,c2 - c21) '
* newc222 = c222 but zero (0) instead of UNDEF values
' define newc222 = const(c2,0,-u)'


* c223 = UNDEF if c212 < 0
' define c223 =  maskout(c21,c21 - c2)'
* newc223 = c223 but zero (0) instead of UNDEF values
' define newc223 = const(c21,0,-u)'


' define d1 = newc122 + newc123 +  newc222 + newc223'
' define d2 = a1 + b1'


' e1 = aave(d1,lon='lon1',lon='lon2m1',lat='lat1',lat='lat2m1')'
' e2 = aave(d2,lon='lon1',lon='lon2m1',lat='lat1',lat='lat2m1')'
' s1 = 100*e2/e1' 

' set gxout print'

' set prnopts %g 1 1'


res = write('outve'r%2004%dir%dt%_%h'.txt','Verif Scores: FCST+'h' hr Init Date '2004%dir%dt)

res = write('outve'r%2004%dir%dt%_%h'.txt','RMSE =')
'd sker'
ret = result
line2 = sublin(ret,2)
number=subwrd(line2,1)
res = write('outve'r%2004%dir%dt%_%h'.txt', number)


res = write('outve'r%2004%dir%dt%_%h'.txt','TEND CORR =')
'd rpr'
ret = result
line2 = sublin(ret,2)
number=subwrd(line2,1)
res = write('outve'r%2004%dir%dt%_%h'.txt', number)


res = write('outve'r%2004%dir%dt%_%h'.txt','MEAN ABS ERR =')
'd  smer'
ret = result
line2 = sublin(ret,2)
number=subwrd(line2,1)
res = write('outve'r%2004%dir%dt%_%h'.txt',number)


res = write('outve'r%2004%dir%dt%_%h'.txt','BIAS =')
'd  bias'
ret = result
line2 = sublin(ret,2)
number=subwrd(line2,1)
res = write('outve'r%2004%dir%dt%_%h'.txt', number)

res = write('outve'r%2004%dir%dt%_%h'.txt','S1 =')
'd  s1'
ret = result
line2 = sublin(ret,2)
number=subwrd(line2,1)
res = write('outve'r%2004%dir%dt%_%h'.txt', number)
'reinit'
h=h+24
endwhile