************************
** Re-write of the plumes script.
* Because I deleted it like an idiot.
************************
'reinit'
say '##############################'
say '##      plot.plumes.gs      ##'
say '##############################'
datefile="date.index"
modfile="grads.models.index"
plumelocfile="grads.locations.plumes.index"
varfile="grads.plumes.index"
outsizex=768
outsizey=1024
say 'Reading 'datefile
rc = 0
while (rc = 0)
   indata=read(datefile)
   rc=sublin(indata,1)
   say rc
   if (rc=0)
    data=sublin(indata,2)
    dtg = subwrd(data,1)
   endif
endwhile
say dtg
*
say 'Reading 'modfile
oo=1
rc = 0
while (rc = 0)
   indata=read(modfile)
   rc=sublin(indata,1)
   say rc
   if (rc = 0)
    data=sublin(indata,2)
    say data
    modtitle.oo = subwrd(data,1)
    label.oo    = subwrd(data,2)
    imtitle.oo  = subwrd(data,3)
    dtype.oo    = subwrd(data,4)
    say oo' 'modtitle.oo' 'label.oo' 'imtitle.oo' 'dtype.oo
   endif
   if (rc = 2)
    nmod = oo-1
   endif
   oo=oo+1
endwhile
say 'nmod 'nmod
*
nn=1
say 'Reading 'plumelocfile
rc = 0
while (rc = 0)
   indata=read(plumelocfile)
   rc=sublin(indata,1)
   say rc
   if (rc = 0)
     data=sublin(indata,2)
     say data
     plat.nn = subwrd(data,1)
     plon.nn = subwrd(data,2)
     say plat.nn' 'plon.nn
   endif
   if (rc=2)
    nloc=nn-1
   endif
   nn = nn + 1
endwhile
*
pl=1
say 'Reading 'varfile
rc=0
while(rc=0)
   indata=read(varfile)
   rc=sublin(indata,1)
   say rc
   if (rc = 0)
     data=sublin(indata,2)
     say data
     varname.pl = subwrd(data,1)
     level.pl   = subwrd(data,2)
     t0.pl      = subwrd(data,3)
     tfin.pl    = subwrd(data,4)
   endif
   if(rc=2)
     npl = pl - 1
   endif
   pl=pl+1
endwhile
say npl
*
c='2 3 4 5 6 7 8 9 10 11 12 13 14'
ccols=''
ii=1
while(ii<=6)
  ccols=ccols' 'c
  ii = ii + 1
endwhile
tt=1
xlabber=''
*tmax=65
*while(tt < tmax/2)
* xlabber=xlabber%' |'
* tt= tt + 1
*endwhile
oo=1
while(oo<=nmod)
  pl=1
  while(pl<=npl)
   tt=1
   xlabber=''
   while(tt < tfin.pl/2)
    xlabber=xlabber%' |'
    tt= tt + 1
   endwhile
   if(level.pl = 'onelev')
    'open onelev.'modtitle.oo'.'dtg'.ctl'
    if(varname.pl = 'wind_10m' | varname.pl = 'wind_19d5m')
     'q file'
     rec = sublin(result,5)
     tmax = subwrd(rec,12)
     zmax = subwrd(rec,9)
     emax = subwrd(rec,15)
     ymax = subwrd(rec,6)
     xmax = subwrd(rec,3)
     'set t 't0.pl' 'tfin.pl
     'set z 1'
     'set y 1 'ymax
     'set x -180 'math_nint(xmax)
     'set e 1 'emax
     if(varname.pl = 'wind_10m')
      'define wind_10m = mag(wnd_ucmp_10m,wnd_vcmp_10m)'
      usevar=varname.pl
     else
      'define wind_19d5m = mag(wnd_ucmp_19d5m,wnd_vcmp_19d5m)'
      usevar=varname.pl
     endif
    endif
    usevar=varname.pl'.1'
   else
    if(varname.pl = 'thickness' )
     'open geop_ht.'modtitle.oo'.'dtg'.ctl'
     'q file'
     rec = sublin(result,5)
     tmax = subwrd(rec,12)
     zmax = subwrd(rec,9)
     emax = subwrd(rec,15)
     ymax = subwrd(rec,6)
     xmax = subwrd(rec,3)
     'set t 't0.pl' 'tfin.pl
     'set z 1'
     'set y 1 'ymax
     'set x -180 'math_nint(xmax/2)
     'set e 1 'emax
     'define thickness = geop_ht(lev=500) - geop_ht(lev=1000)'
     usevar=varname.pl
    else
     if(varname.pl = 'wind')
      'open wnd_ucmp.'modtitle.oo'.'dtg'.ctl'
      'q file'
      rec = sublin(result,5)
      tmax = subwrd(rec,12)
      zmax = subwrd(rec,9)
      emax = subwrd(rec,15)
      ymax = subwrd(rec,6)
      xmax = subwrd(rec,3)
      'set t 't0.pl' 'tfin.pl
      'set lev 'level.pl
      'set y 1 'ymax
      'set x 1 'xmax
      'set e 1 'emax
      'open wnd_vcmp.'modtitle.oo'.'dtg'.ctl'
      'define wind = mag(wnd_ucmp.1,wnd_vcmp.2)'
      'close 2'
      usevar="wnd_ucmp"
     else
      'open 'varname.pl'.'modtitle.oo'.'dtg'.ctl'
      usevar=varname.pl'.1'
      'q file'
      rec = sublin(result,5)
      tmax = subwrd(rec,12)
      zmax = subwrd(rec,9)
      emax = subwrd(rec,15)
      ymax = subwrd(rec,6)
      xmax = subwrd(rec,3)
      'set t 't0.pl' 'tfin.pl
      'set lev 'level.pl
      'set y 1 'ymax
      'set x 1 'xmax
      'set e 1 'emax
     endif
    endif
   endif
   if(level.pl = 'onelev')
    'open ensmean.onelev.'varname.pl'.'modtitle.oo'.'dtg'.ctl'
   else
    'open ensmean.'varname.pl'.'modtitle.oo'.'dtg'.ctl'
   endif
   'set t 1'
   'set e 1'
   'set y 1 'ymax
   'set x 1 'tfin.pl
   'define median = 'usevar
   'define pc10 = 'usevar
   'define pc90 = 'usevar
   'define pc25 = 'usevar
   'define pc75 = 'usevar
*
   nn=1
   'q ctlinfo 1'
   rec=sublin(result,3)
   undef=subwrd(rec,2)
   while(nn<=nloc)
    'set lat 'plat.nn
    'set lon 'plon.nn
    'set e 1 last'
    'set t 1 'tfin.pl
    if (level.pl = 'onelev')
     'set z 1'
    else
     if (varname.pl = 'thickness')
      'set z 1'
     else
      'set lev 'level.pl
     endif
    endif
    'set gxout stat'
    'd 'varname.pl
    rec = sublin(result,9)
    vrmin=subwrd(rec,5) - subwrd(rec,7)
    vrmax=subwrd(rec,6) + subwrd(rec,7)
    vrdiff=vrmax-vrmin
    say vrdiff
    if(vrdiff >= 1000)
     extra=100
    else
     if(vrdiff >= 500)
      extra=50
     else
      if(vrdiff >= 300)
       extra=25
      else
       if(vrdiff >= 100)
        extra=15
       else
        if(vrdiff >= 50)
         extra=5
        else
         if(vrdiff >= 10)
          extra=2.5
         else
          if(vrdiff >= 5)
           extra=0.5
          else
           if(vrdiff >= 1)
            extra=0.25
           else
            if(vrdiff >= 0.1)
             extra=0.05
            else
             extra=0.01
            endif
           endif
          endif
         endif
        endif
       endif
      endif
     endif
    endif
    if(varname.pl = 'ttl_prcp_06' | varname.pl = 'ttl_prcp')
     vrmin=0
    endif
    'set vpage 0 8.5 5.5 11'
    'set parea 1 7.5 1.00 4.75'
    'set grads off'
    'set gxout contour'
    'set vrange 'vrmin' 'vrmax
    'set ylint 'extra
    'set cmark 0'
    ee=1
    while(ee<=emax)
     'set e 'ee
     cc=subwrd(ccols,ee)
     'set ccolor 'cc
     'set cmark 0'
     if(ee=emax)
      'set xlabs'xlabber
     endif
     'set xlab on'
     'd 'varname.pl
     ee = ee + 1
    endwhile
    'q file 2'
    line=sublin(result,7)
    units=subwrd(line,4)
    'plumetitle.gs 'label.oo' 'dtg' 'varname.pl' 'level.pl' 'units' 'plat.nn' 'plon.nn
    'draw xlab Valid Time (UTC)'
    'set string 1'
    'set strsiz 0.11 0.12'
    'draw string 7.6 4.65 Members'
    'set line 1 1'
    'draw line 7.55 4.6 8.45 4.6'
**
* Getting number of member models.
**
    'q file 2'
    line=sublin(result,6)
    nvars=subwrd(line,5)
    m=1
    while(m<=nvars)
     ii = 6 + m
     'q file 2'
     line=sublin(result,ii)
     var=subwrd(line,1)
     if(var=sd)
      jj=ii
      break
     endif
     m=m+1
    endwhile
    nens=nvars-m
    n=1
    ystart=4.4
    yincr=0.2
    'set e 1'
    while (n<=nens)
     ii = jj + n
     'q file 2'
     line=sublin(result,ii)
     say line
     var=subwrd(line,1)
     ensname=subwrd(line,4)
     'q dims 2'
     say result
     'set gxout stat'
     say varname.pl
     'd 'var'.2'
     say result
     line=sublin(result,8)
     ensn=subwrd(line,5)
     'set strsiz 0.11 0.12'
     ypos = ystart - yincr*(n-1)
     'draw string 7.6 'ypos' 'ensname
     'draw string 8.1 'ypos' 'ensn
     n=n+1
    endwhile
    'set gxout contour'
    tt=1
    while(tt<=tfin.pl)
     emaxtrue=emax
     ee=1
     'set t 'tt
     while(ee<=emax)
      'set e 'ee
      'd 'usevar
      line=sublin(result,1)
      val.ee=subwrd(line,4)
      if(val.ee=undef)
       emaxtrue=emaxtrue-1
      endif
      ee = ee + 1
     endwhile
     ee=1
     while(ee<=emax)
      temp=-999999
      ff=emax
      while(ff>ee)
       gg=ff-1
       if(val.gg < val.ff)
        temp = val.gg
        val.gg = val.ff
        val.ff = temp
       endif
       ff = ff - 1
      endwhile
      if (temp = -999999) ; break ; endif
      ee = ee + 1
     endwhile
     lower=math_int(0.9*emaxtrue)
     upper=math_int(0.1*emaxtrue)+1
     prc90.tt=val.upper
     prc10.tt=val.lower
     lower=math_int(0.75*emaxtrue)
     upper=math_int(0.25*emaxtrue)+1
     prc75.tt=val.upper
     prc25.tt=val.lower
     midcheck=math_mod(emaxtrue,2)
     if(midcheck = 0)
      upper=math_int(0.5*emaxtrue)
      lower=math_int(0.5*emaxtrue)+1
      median.tt=(val.upper+val.lower)/2
     else
      upper=math_nint(0.5*emaxtrue)
      median.tt=val.upper
     endif
*     if(tt=1)
*      ee=1
*      while (ee<=emax)
*       say val.ee
*       ee=ee+1
*      endwhile
*      say
*      say prc90.tt
*      say prc10.tt
*     endif
     tt = tt + 1
    endwhile
    'q dims 1'
    rec=sublin(result,3)
    yval=subwrd(rec,9)
    say yval
    tt = 1
    while(tt<=tfin.pl)
     say tt' 'yval' 'prc90.tt' 'prc75.tt' 'median.tt' 'prc25.tt' 'prc10.tt
     'set defval pc90 'tt' 'yval' 'prc90.tt
     'set defval pc75 'tt' 'yval' 'prc75.tt
     'set defval median 'tt' 'yval' 'median.tt
     'set defval pc25 'tt' 'yval' 'prc25.tt
     'set defval pc10 'tt' 'yval' 'prc10.tt
     tt = tt + 1
    endwhile
    'set lat 'plat.nn
    'set lon 'plon.nn
    'set e 1'
    'set t 1 'tfin.pl
    'define ensmin=tloop(min('usevar',e=1,e='emax'))'
    'define ensmax=tloop(max('usevar',e=1,e='emax'))'
    'set vpage 0 8.5 2.7 6'
    'set parea 1 7.5 0.7 2.8'
    'set grads off'
    'set t 1 'tfin.pl
    'set vrange 'vrmin' 'vrmax
    'set gxout linefill'
    'set lfcols 1 5'
    'set xlpos 0 t'
    'set ylint 'extra*2
    'd ensmin ; ensmax'
    'set xlab off'
    'set t 1'
    'set x 1 'tfin.pl
    'set gxout line'
    'set xlab off'
    'set cmark 0'
    'set ccolor 1'
    'set cthick 4'
    'set cstyle 1'
    'd pc90'
    'set xlab off'
    'set cmark 0'
    'set ccolor 1'
    'set cthick 4'
    'set cstyle 1'
    'd pc10'
    'set xlab on'
    'set gxout linefill'
    'set lfcols 0 7'
    'set xlabs'xlabber
    'd pc25 ; pc75'
    'set gxout line'
    'set xlab off'
    'set cmark 0'
    'set ccolor 2'
    'set cthick 6'
    'set cstyle 1'
    'd median'
*
    'set strsiz 0.11 0.12'
    'set string 2'
    'draw string 7.6 2.7 Median'
    'set string 7'
    'draw string 7.6 2.5 75%-25%'
    'set string 1'
    'draw string 7.6 2.3 90%-10%'
    'set string 5'
    'draw string 7.6 2.1 Max/Min'
*
    'set vpage 0 8.5 0 3.3'
    'set parea 1 7.5 1 3.1'
    'set grads off'
    'set lat 'plat.nn
    'set lon 'plon.nn
    'set vrange 'vrmin' 'vrmax
    'set ylint 'extra*2
    'set xlab on'
    'set t 1 'tfin.pl
    'set gxout errbar'
    'set ccolor 4'
    'set xlpos 0 b'
    'd ensmin ; ensmax'
    'define plus = mean.2 + sd.2'
    'define minus = mean.2 - sd.2'
    'set gxout bar'
    'set bargap 50'
    'set baropts outline'
    'set ccolor 3'
    'd minus ; plus'
    'set gxout line'
    'set xlabs'xlabber
    'set cmark 0'
    'set ccolor 6'
    'set cthick 6'
    'set cstyle 1'
    'd mean.2'
    'draw xlab Valid Time (UTC)'
*
    'set strsiz 0.11 0.12'
    'set string 6'
    'draw string 7.6 2.99 Mean'
    'set string 3'
    'draw string 7.6 2.79 Std. Dev.'
    'set string 4'
    'draw string 7.6 2.59 Max/Min'
*
    if(plat.nn = 0)
     uselat='0'
    else
     if(plat.nn > 0)
      uselat=plat.nn%'N'
     else
      uselat=(-1)*plat.nn%'S'
     endif
    endif
    if(plon.nn = 0)
     uselon='0'
    else
     if(plon.nn > 0)
      uselon=plon.nn%'E'
     else
      uselon=(-1)*plon.nn%'W'
     endif
    endif
    if(level.pl = 'onelev')
     'gxyat -x 'outsizex' -y 'outsizey' 'imtitle.oo'.plume-'varname.pl'.'uselat'.'uselon'.'dtg'.png'
    else
     if(varname.pl = 'thickness')
      'gxyat -x 'outsizex' -y 'outsizey' 'imtitle.oo'.plume-'varname.pl'.'uselat'.'uselon'.'dtg'.png'
     else
      'gxyat -x 'outsizex' -y 'outsizey' 'imtitle.oo'.plume-'varname.pl'_'level.pl'mb.'uselat'.'uselon'.'dtg'.png'
     endif
    endif
    'c'
    nn = nn + 1
   endwhile
   pl=pl+1
   'close 2'
   'close 1'
  endwhile
  oo=oo+1
endwhile