function exce_ver (argument_list)
*
* count where excedence is showing the right sign
* comparing integrated exceedance duration with
* integrated (4weeks) Precip, T, ff10.

* time series at one location (may be area mean) for one variable
*
* input: obs.ctl   connecting to ecmwf 12-36 hrs forecast data
*        mclemean.ctl hindcast average from several forecasts in formost archive
*        mclmedia.ctl hindcast median  from several forecasts in formost archive
*        lat
*        lon
*        parameter= one of 'tmean tmax tmin prec sun mslp ff10'
*
* output: plots of fraction above / below forecasts
*
'reinit'
'run /home/fr2100/frtm/bin/grads/becol.gs'

'!echo $base'
base=result              

say base

say '  call with excedence_ver.gs fxdate area var'
say ' '
say 'allowed parameters for area: uk europe global '
say 'allowed parameters for  var: heat drought flood cold winds      '

base='/data/nwp1/frtm/formost'
base='/project/formost/'
obase='/data/local/frtm/bormost'

say base
*fxdate='20040421'

fxdate=subwrd(argument_list,1)
area  =subwrd(argument_list,2)
var   =subwrd(argument_list,3)
say argument_list


* calculate time window with respect to fxdate
*  we have "heat exceedance products" since 20070510
*  we have "flood and drought exceedance products" since 20070531
*  we have "cold  and wind    exceedance products" since 20080501

tag=var
if ( var = heat )                  
     startdate = 20070510 
     startdate = 20070614 
     tag=''
     par=tmax
*   to compare with monthly means Obs!
     thresh=22
endif
if ( var = flood ) 
     startdate = 20070531 
     startdate = 20070621 
     par=prec
     thresh=3 
endif
if ( var = drought ) 
     startdate = 20070614 
     par=prec
     thresh= 1
endif
if ( var = winds ) 
     startdate = 20080501 
     par=ff10
     thresh=5
endif
if ( var = cold    ) 
     startdate = 20071025 
     par=tmax
     thresh=5
endif

if ( area = 'global' | area = '' )

   lat1=-40
   lat2=70
   lon1=-130
   lon2=160

endif

if ( area = 'europe' )

   lat1=30
   lat2=60
   lon1=-20
   lon2=40

endif

if ( area = 'uk' )

   lat1=49
   lat2=60
   lon1=-12
   lon2=4

endif

say 'datediff  'fxdate' 'startdate' >./erp2'
'!datediff  'fxdate' 'startdate' >./erp2'

record2=read  (erp2)
timenode=sublin(record2,2)
rc     =close (erp2)
maxfx = timenode/7
*maxfx = 26             

outdir=base'/../public_html/formost/verification/'
outdir='/home/h02/frtm/public_html/formost/verification/'
 
*say outdir 

'set gxout grid'

* loop over the forecasts

fxnum=0

while ( fxnum < maxfx )
   fxnum=fxnum +1

   delta=  (fxnum -1)*7
   delta='-- -'delta
   bxdate=newdate(fxdate,delta)

* 1 89 - 98 era40 mean
       'open 'base'/'bxdate'/mobsemean.ctl'
*say '1: open 'base'/'bxdate'/mobsemean.ctl'

* 2 90 member 18 year hindcast
       'open 'base'/'bxdate'/mclmedia.ctl'
*say '2: open 'base'/'bxdate'/mclmedia.ctl'

* 3    verification data from 12 - 36 hrs high resolution forecast 
       'open 'base'/'bxdate'/mobs.ctl'
*say '3: open 'base'/'bxdate'/mobs.ctl'

* 4   ensemble mean forecast
       'open 'base'/'bxdate'/mfxemean.ctl'
*say '4: open 'base'/'bxdate'/mfxemean.ctl'

* 5
       'open 'base'/'bxdate'/'tag'pdffx.ctl'
*say '5: open 'base'/'bxdate'/'tag'pdffx.ctl'

* 6
       'open 'base'/'bxdate'/'tag'pdfhc.ctl'
*say '6: open 'base'/'bxdate'/'tag'pdfhc.ctl'

* 7 89 - 98 sigma
       'open 'base'/'bxdate'/mobssigma.ctl'
*say '7: open 'base'/'bxdate'/mobssigma.ctl'

* 8 hindcast
       'open 'base'/'bxdate'/mclemean.ctl'
*say '8: open 'base'/'bxdate'/mclemean.ctl'

   'set grads off'
   'set datawarn off'

*'set mproj nps'
   'set lon 'lon1' 'lon2
   'set lat 'lat1' 'lat2
*   say "lat lon " lat1' 'lat2' 'lon1' 'lon2

   'set dfile 1'
   'set t 1'

* monthly mean observed 90ties climatology

   'define mmobsemean  = ('par'.1(t=4)+'par'.1(t=11)+2.*'par'.1(t=18))/4'
   'define thresh=const(mmobsemean,'thresh',-a)'
   'define mmobsigma= ('par'.7(t=4)+'par'.7(t=11)+2.*'par'.7(t=18))/4'

* observed monthly mean

   'define mmobs  =  'par'.3(t=2)'

* forecast monthly ensemble mean

   'define mmobsfx  =  'par'.4(t=2)'

* hindcast monthly ensemble mean

   'define mmobscl  =  'par'.8(t=2)'

* root mean square distance between er40 89- 98 and model hindcast             

   'define rmsdi  =   pow(mmobscl - mmobsemean,2)       '

* root mean square distance between observed monthly mean and forecast ensemble mean
*  'define rmsdi  =   pow(mmobsfx - mmobs,2)       '

   'define rmsdi  = const(rmsdi,0,-u) '

* predicted days above threshold

   'define predex=const('par'.5(t=3),0,-a)'
*  'd predex'

   tim=2
   while ( tim < 32 )

      'define predex = predex + 'par'.5(t='tim')/100. '
* say 'define predex = predex + 'par'.5(t='tim')/100. '
      tim=tim+1
   endwhile

* integrated climatology of   days above threshold

   'define climex=const('par'.6(t=3),0,-a)'

* ascertain that the first time slice is defined, else all climex will be undefined. 
   tim=2
   while ( tim < 32 )

      'define climex = climex + 'par'.6(t='tim')/100. '
      tim=tim+1
   endwhile

* check where the observed monthly mean was larger than the
* climate monthly mean by at least one/half standard deviation, true=1, false=0

*  'define  mmobs=         mmobs+0.5 * mmobsigma'
*  'define  mmobs=         mmobs+thresh'
   'define  event= maskout(mmobs,thresh-mmobs)'
   'define    v1=const(const(event,0),1,-u)'
   'define notv1=const(const(event,1),0,-u)'

    rel1 = '>'
    rel2 = '+'
    if ( var = cold | var = drought ) 
*                    invert test for cold, i.e. associate true with observed < climatology
           'define    v1 = pow(   v1 - 1,2)'
           'define notv1 = pow(notv1 - 1,2)'
            rel1 = '<'
            rel2 = '-'
    endif

* check where the integrated predicted exceedances are larger than the
* climate exceedances , true=1, false=0
    percentage = 30

*   say 'define event=maskout(predex,predex+0.'percentage'-climex)'
   'define event=maskout(predex,predex+0.'percentage'-climex)'
   'define    v2=  const(const( event ,0),1,-u)'
   'define notv2=  const(const( event ,1),0,-u)'

* check where the predicted threshold exceedance is observed. 

    if ( fxnum = 1 )
      'define cnt=       v1'
      'define hit=       v1*v2'                      
      'define fa =    notv1*v2'                     
      'define mis=       v1*notv2'                   
      'define cor=    notv1*notv2'                  
      'define rmsds  =       rmsdi'                   
    else
      'define cnt=   cnt +    v1'                
      'define hit=   hit +    v1*v2'                
      'define fa =   fa  + notv1*v2'                
      'define mis=   mis +    v1*notv2'             
      'define cor=   cor + notv1*notv2'             
      'define rmsds= rmsds + rmsdi'                   
* integrated predicted days above threshold
    endif

*pull a
    'close 8'
    'close 7'
    'close 6'
    'close 5'
    'close 4'
    'close 3'
    'close 2'
    'close 1'

endwhile  * from while over forecasts fxnum

say maxfx'  'fxnum
'define a=hit'
'define b=fa'
'define c=mis'
'define d=cor'
say 'define fxnum=const(a,'maxfx',-a)'
'define fxnum=const(a,'maxfx',-a)'
say 'define s= (a+c)/fxnum'
'define s= (a+c)/fxnum'
say 'define H  = a /( a+c)'
'define H  = const (a /( a+c)     ,0,-u)'
say 'define F  = b /( b+d)'
'define F  = const (  b /( b+d)     ,0,-u)'
say 'define hss= (2*s*(1-s)*(H-F))/(s+s*(1-2*s)*H+(1-s)*(1-2*s)*F)'
'define hss= (2*s*(1-s)*(H-F))/(s+s*(1-2*s)*H+(1-s)*(1-2*s)*F)'
'define hss= const ( hss , 0, -u)'
say 'define pss= (a*d-b*c)/((b+d)*(a+c))'
'define pss= const ( (a*d-b*c)/((b+d)*(a+c)) , 0, -u)'
'define csi= const ( a/(b+a+c) , 0, -u)'
'define YQ = const ( (H-F)/ (H*(1-F)+F*(1-H)) ,0,-u)'
'define OR = const ( H/(1-H) * (1-F)/F        ,0,-u)'
 
say 'define rmsd= sqrt(rmsds/'maxfx')'
'define rmsd= sqrt(rmsds/'maxfx')'

say ' define EDS=2*log((a+c)/fxnum)/log(a/fxnum) - 1.'
'define EDS=       2*log((a+c)/fxnum)/log(a/fxnum) - 1.'
'define EDS=const( EDS ,0,-u)'

*pull x

* second part of header string:
part2='\ Observed value 'rel1' 90s climate 'rel2' 1/2 sigma and \ Area under the fx PDF > than Area under the hindcast PDF'
part2='\ Observed monthly mean 'rel1' 'thresh ' and \ Area under the fx PDF 'percentage'% > than Area under the hindcast PDF'

'c'
'set grads off'
*'set gxout grid'
'set gxout contour'
'd cnt'
'draw map'
'draw title observed count 'var' 'bxdate' to 'fxdate'  'part2
printname='thver_'fxdate var area'_cnt' 
say 'printim 'outdir'/'printname'.gif gif white'
'printim 'outdir'/'printname'.gif gif white'

'c'
'set grads off'
'set gxout contour'
*'set gxout shaded'
'd pss'
'draw map'
'draw title Peirce skill score 'var' 'bxdate' to 'fxdate'  'part2
printname='thver_'fxdate var area'_peirce' 
say 'printim 'outdir'/'printname'.gif gif white'
'printim 'outdir'/'printname'.gif gif white'

'c'
'set grads off'
'set gxout contour'
*'set gxout shaded'
'd csi'
'set clevs 0.5'
'set ccols   1 '
'set cthick 6'
'd csi'
'draw map'
'draw title Critical success rate 'var' 'bxdate' to 'fxdate'  'part2 
printname='thver_'fxdate var area'_csi' 
say 'printim 'outdir'/'printname'.gif gif white'
'printim 'outdir'/'printname'.gif gif white'

'c'
'set grads off'
'set gxout contour'
*'set gxout shaded'
'd  YQ'
'draw map'
'draw title Yules Q 'var' 'bxdate' to 'fxdate'  'part2
printname='thver_'fxdate var area'_yq' 
say 'printim 'outdir'/'printname'.gif gif white'
'printim 'outdir'/'printname'.gif gif white'

'c'
'set grads off'
'set gxout contour'
*'set gxout shaded'
'd hss'
'draw map'
'draw title Heidke skill score 'var' 'bxdate' to 'fxdate' 'part2
printname='thver_'fxdate var area'_heidke' 
say 'printim 'outdir'/'printname'.gif gif white'
'printim 'outdir'/'printname'.gif gif white'

'c'
'set grads off'
'set gxout contour'
*'set gxout shaded'
'd rmsd'
'draw map'
'draw title RMS Distance  'par' 'bxdate' to 'fxdate' between \ hindcast and ERA40 89 to 98 '

printname='thver_'fxdate var area'_rmsd' 
say 'printim 'outdir'/'printname'.gif gif white'
'printim 'outdir'/'printname'.gif gif white'

quit

  function newdate(arg1,arg2,arg3,arg4)
*
*usage: newdate [-s] yyyymmddhh [+|-]hours
*       newdate [-s] -D yyyymmdd [+|-]days
*       newdate [-d | -j | -J | -s] yyyymmddhh
*       -D  switch to use days instead of hours
*       -d  prints the Julian date
*       -j  prints the julian day number of the year
*       -J  prints the julian day number of the year (JJJ)
*       -s  prints the yyyymmhhdd separately (yyyy mm hh dd)
*       Julian number is in range of [1-366]
*
  '!newdate -D 'arg1' 'arg2' 'arg3' 'arg4' >./erp1'
*  say '!newdate -D 'arg1' 'arg2' 'arg3' 'arg4
  record2=read(erp1)
  result=sublin(record2,2)
  rc= close (erp1)
  return(result)