'reinit'
'open emc_statistics.ctl'
'set xlopts 1 9 0.20'
'set ylopts 1 9 0.20'
*'set parea 1.8 9.5 1.8 7'
'set parea 1.0 10.2 1.0 7'
*
* map the model with the fort.11 file
*
rc=read('fort.11')
ptype=sublin(rc,2)
fname='st05'
say "Type of plot is "fname
rc=read('fort.11')
ntime=sublin(rc,2)
say "Number of verification time is "ntime
rc=read('fort.11')
nmodel=sublin(rc,2)
say "Number of models to process is "nmodel
rc=read('fort.18')
legend_opt=sublin(rc,2)
if(legend_opt='');legend_opt=1;endif
rc=read('fort.12')
exp_tilte=sublin(rc,2)
rc=read('fort.12')
prefix=sublin(rc,2)
exp_flag=substr(prefix,3,2)
rc=read('fort.12')
iline=sublin(rc,2)
errbar=subwrd(iline,2)
imodel=1
while (imodel <= nmodel)
 rc=read('fort.12')
 eof=sublin(rc,1)
 if (eof=0)
  iline=sublin(rc,2) 
  model.imodel=subwrd(iline,1)
  cmark.imodel=subwrd(iline,2)
  ctype.imodel=subwrd(iline,3)
  ccolo.imodel=subwrd(iline,4)
  cthik.imodel=subwrd(iline,5)
 else
  say "error reading file...too many models "eof
  break
 endif
 rc1=read('fort.14')
 eof=sublin(rc1,1)
 if (eof=0)
  note.imodel=sublin(rc1,2)
 else
  say "error reading file...too many models "eof
  break
 endif
 say "Model to prosses is "model.imodel
 say "cmark to prosses is "cmark.imodel
 say "ctype to prosses is "ctype.imodel
 say "ccolo to prosses is "ccolo.imodel
 say "cthik to prosses is "cthik.imodel
 say "Exp flag is "exp_flag
 say "Note for model is "note.imodel
 imodel = imodel+1
endwhile
*
* set the mark for the OFCL time
*J.Peng------
nt.1=1
nt.2=3
nt.3=5
nt.4=7
nt.5=9
nt.6=11
nt.7=13
nt.8=15
nt.9=17
nt.10=19
nt.11=21
nt.12=23
nt.13=25
nt.14=27
nt.15=29

*
* plotting error
*
'set strsiz 0.14'
'set grads off'
'set lev 1'
*J.Peng 'set xaxis 0 120 12'
'set xaxis 0 168 12'

'set ylint 5'
'set gxout line'
*'set string 1 c 9'
*'draw string 10.4 6.94 LEGEND:'
ilev=1
cmaxv=25
cminv=-25
dint=(cmaxv-cminv)/5
while(ilev<=nmodel)
 'set vrange 'cminv' 'cmaxv
 'set lev 'ilev
 'set ccolor 'ccolo.ilev
 'set cthick 'cthik.ilev
 'set cmark 0'
 'set cstyle 'ctype.ilev
 'set ylint 'dint
 if (fname="m34")
  'define m34=(ne34+se34+nw34+sw34)/4'
 endif
 if (fname="m50")
  'define m50=(ne50+se50+nw50+sw50)/4'
 endif
 if (fname="m65")
  'define m65=(ne65+se65+nw65+sw65)/4' 
 endif
 if (fname="st01")
  'define st01=-(tker-ave(tker,lev=1,lev=1))*100/ave(tker,lev=1,lev=1)'
  relative_model=model.1
 endif
 if (fname="st02")
  'define st02=-(tker-ave(tker,lev=2,lev=2))*100/ave(tker,lev=2,lev=2)'
  relative_model=model.2
 endif
 if (fname="st03")
  'define st03=-(tker-ave(tker,lev=3,lev=3))*100/ave(tker,lev=3,lev=3)'
  relative_model=model.3
 endif
 if (fname="st04")
  'define st04=-(tker-ave(tker,lev=4,lev=4))*100/ave(tker,lev=4,lev=4)'
  relative_model=model.4
 endif
 if (fname="st05")
  'define st05=-(tker-ave(tker,lev=5,lev=5))*100/ave(tker,lev=5,lev=5)'
  relative_model=model.5
 endif
 if (fname="si01")
  'define si01=-(wind-ave(wind,lev=1,lev=1))*100/ave(wind,lev=1,lev=1)'
  relative_model=model.1
 endif
 if (fname="si02")
  'define si02=-(wind-ave(wind,lev=2,lev=2))*100/ave(wind,lev=2,lev=2)'
  relative_model=model.2
 endif
 if (fname="si03")
  'define si03=-(wind-ave(wind,lev=3,lev=3))*100/ave(wind,lev=3,lev=3)'
  relative_model=model.3
 endif
 if (fname="si04")
  'define si04=-(wind-ave(wind,lev=4,lev=4))*100/ave(wind,lev=4,lev=4)'
  relative_model=model.4
 endif
 if (fname="si05")
  'define si05=-(wind-ave(wind,lev=5,lev=5))*100/ave(wind,lev=5,lev=5)'
  relative_model=model.5
 endif
 'd 'fname
*
* draw mark
*
 icount=1
*J.Peng  while (icount <=12)
while (icount <=15)

  ntim=nt.icount
  'd ave('fname',x='ntim',x='ntim')'
  rc=sublin(result,2)
  value=subwrd(rc,4)
  'q gr2xy 'ntim' 'value
  rc=sublin(result,1)
  xc=subwrd(rc,3)
  yc=subwrd(rc,6)
  'set line 'ccolo.ilev
  'draw mark 'cmark.ilev' 'xc' 'yc' 0.16'
*
* reserved the capability of drawing bias/radii error bars for testing now. Need to pull
* error bar information for these variables from the verify_model.f. 
*
  if (errbar="Y" | errbar="y")
   if (fname='wind' |fname='biass' |  fname='tker' | fname='m34' | fname='m50' | fname='m65')
    if (fname="bias")
     'd ave(wind_std,x='ntim',x='ntim')'
    else
     'd ave('fname'_std,x='ntim',x='ntim')'
    endif
    rc=sublin(result,2)
    std=subwrd(rc,4)
    valuep=value + 1.96*std
    valuem=value - 1.96*std
    'q gr2xy 'ntim' 'valuep
    rc=sublin(result,1)
*J.Peng    xp=subwrd(rc,3)+(ilev-1)*0.00
    xp=subwrd(rc,3)+(ilev-1)*0.005

    yp=subwrd(rc,6)
    'q gr2xy 'ntim' 'valuem
    rc=sublin(result,1)
*J.Peng    xm=subwrd(rc,3)+(ilev-1)*0.00
    xm=subwrd(rc,3)+(ilev-1)*0.005

    ym=subwrd(rc,6)
    'set line 'ccolo.ilev' 'ctype.ilev' 'cthik.ilev
    'draw line 'xm' 'ym' 'xp' 'yp
    'draw line 'xm-0.05' 'ym' 'xm+0.05' 'ym
    'draw line 'xp-0.05' 'yp' 'xp+0.05' 'yp
   endif
  endif
  icount=icount+1
 endwhile
 ilev=ilev+1  
endwhile
*
* draw labels
*
ix=2.0
iy=7.2
iy=iy+0.63
'set string 1 c 9 90'
'draw string 0.36 4.0 TRACK FORECAST SKILL (%)'
'set string 1 c 9 0'
'draw string 5.5 0.25 Forecast lead time (hr)'
'set string 1 c 9 0'
'set strsiz 0.18'
if (exp_flag = 'AL' |  exp_flag = 'al') 
 'draw string 5.8 'iy-0.3' 'exp_tilte
else
 'draw string 5.8 'iy-0.3' STATISTISCS FOR A SINGLE CASE - 'prefix
endif
*J.Peng 'draw string  5.8 'iy' HWRF FORECAST - TRACK FORECAST SKILL (%) STATISTICS '
'draw string  5.8 'iy' MODEL FORECAST - TRACK FORECAST SKILL (%) STATISTICS '

if (fname="st01" | fname="st02" | fname="st03" | fname="st04" | fname="st05")
 'set string 1 l 18 0'
 'draw string 1.1 1.3 SKILL PLOT RELATIVE TO THE 'relative_model' MODEL'
endif
if (fname="si01" | fname="si02" | fname="si03" | fname="si04" | fname="si05")
 'set string 1 l 18 0'
 'draw string 1.1 1.3 SKILL PLOT RELATIVE TO THE 'relative_model' MODEL'
endif
*
* add ncase finally
*
'set cmark 0'
'set cthick 5'
'set ccolor 1'
'set lat 0'
'd lat'
'set gxout print'
*J.Peng 'set prnopts %4.0f 21 1'
'set prnopts %4.0f 29 1'

'd nc_track'
record=sublin (result,2)
'set string 5 c 9'
'set strsiz 0.14'
'draw string 0.4 0.5 #CASE'
itime=1
xcase=1.0
ycase=0.5
*J.Peng while (itime <=21)
while (itime <=29)

 ncase=subwrd(record,itime)
 if (ncase > 0)
  'set string 5 c 9'
  'draw string 'xcase' 'ycase' 'ncase
 else
  'set string 0 c 9'
  'draw string 'xcase' 'ycase' . '
 endif
*J.Peng  xcase=xcase+0.92
 xcase=xcase+0.65

 itime=itime+2
endwhile
#
# draw legend
#
if (legend_opt=1)
 ix=2.0
 iy=7.2
else
 ix=2.8
 iy=6.8
endif
ilev=1
while (ilev<=nmodel)
 if (legend_opt=1)
  if (ix > 9.5)
   iy = iy + 0.23
   ix = 2.0
  endif
  'set strsiz 0.11'
  'set string 'ccolo.ilev' l 'cthik.ilev
  'draw string 'ix-0.4' 'iy' 'model.ilev': 'note.ilev
  is=ix-1.0
  ie=ix-0.5
  ic=ix-0.75
  'set line 'ccolo.ilev' 'ctype.ilev' 'cthik.ilev
  'draw line 'is' 'iy' 'ie' 'iy
  'draw mark 'cmark.ilev' 'ic' 'iy' 0.16'
  ix=ix+2.8
 else
  'set strsiz 0.18'
  'set string 'ccolo.ilev' l 20'
  'draw string 'ix' 'iy' 'model.ilev': 'note.ilev
  is=ix-1.6
  ie=ix-0.18
  ic=ix-0.9
  'set line 'ccolo.ilev' 'ctype.ilev' 18'
  'draw line 'is' 'iy' 'ie' 'iy
  'draw line 'is' 'iy-0.04' 'ie' 'iy-0.04
  'draw mark 'cmark.ilev' 'ic' 'iy-0.02' 0.16'
  iy=iy-0.3
 endif
 ilev=ilev+1
endwhile
#
# print title and convert to png file
#
'set string 1 r 5'
'set strsiz 0.1'
*'draw string 10.6 0.2 HWRF project - NOAA/NCEP/EMC'
'enable print fig_emc_'fname'.gmf'
'print'
'disable print'
'!gxps -c -i fig_emc_'fname'.gmf -o fig_emc_'fname'.ps'
'!convert -trim -rotate 90 -geometry 700x700 fig_emc_'fname'.ps fig_'prefix'_'fname'.png'
*
* draw legend
*
*ix=5
*iy=7
*ilev=1
*while (ilev<=nmodel)
* 'set strsiz 0.36'
* 'set string 'ccolo.ilev' l 18'
* 'draw string 'ix' 'iy' 'model.ilev': 'note.ilev
* is=ix-2.5
* ie=ix-0.5
* ic=ix-1.5
* 'set line 'ccolo.ilev' 'ctype.ilev' 18'
* 'draw line 'is' 'iy' 'ie' 'iy
* 'draw line 'is' 'iy-0.04' 'ie' 'iy-0.04
* 'draw mark 'cmark.ilev' 'ic' 'iy-0.02' 0.16'
* iy=iy-0.7
* ilev=ilev+1
*endwhile
*'enable print fig_emc_'fname'.gmf'
*'print'
*'disable print'
*'!gxps -c -i fig_emc_'fname'.gmf -o fig_emc_'fname'.ps'
*'!convert -trim -rotate 90 -geometry 700x700 fig_emc_'fname'.ps fig_'prefix'_legend.png'
'!rm -f *.gmf *.ps'