*
*  Calculates and plots 5-month averages of precipitation
*
*  To use, grads -p
*
'reinit'
*
*  allow GrADS script functions...
*
rc = gsfallow("on")
*
'open pptn.ctl'
'enable print pptn_pentad.gx'
'set lon 0 360'
'set lat -90 90'
'set z 1'
*
*  ts = starting time for averages and plots
*  tf = final time for averages and plots
*
ts = 3 ; tf = 118

*---------------------------------------------
*  Method 1:  define and plot 5-month averages
*---------------------------------------------

'set t 'ts' 'tf
'define pptn5a = ave(pptn,t-2,t+2)'

tt = ts ; n = 1
while ( tt <= tf )
  say
  say 'Plot number 'n
  say '...time levels 'tt-2' 'tt+2
  'c'
  'set t 'tt
  'q time'
  date = subwrd(result,3) ; month = substr(date,6,3)
*
*  define contour levels
*
  'set gxout stat'
  'd pptn'
  min_max = sublin(result,8)
  rec14 = sublin(result,14)
  min = subwrd(min_max,4) ; max = subwrd(min_max,5)
  nmin = math_nint(min+0.5) ; nmax = math_nint(max-0.5)
*  rmin = min+0.5 ; rmax = max-0.5
  say ' min = 'min
  say ' max = 'max
  say ' nmin = 'nmin
  say ' nmax = 'nmax
*  say ' rmin = 'rmin
*  say ' rmax = 'rmax
  pmin = nmin
  pmax = nmax
  ncols = 10 ; nlevs = ncols+1
  int = (pmax-pmin)/ncols ; say 'int = 'int
  clevels = ''
  l=1
  clev = pmin
  while ( l <= nlevs )
    clevels = ''clevels' 'clev
    clev = clev + int
    l=l+1
  endwhile
  say 'clevels = 'clevels
  'set vpage  0 8.5 5.5 11'
  'set parea  1 7   1 4.5'
  'set grads off'
  'set gxout shaded'
  'set clevs 'clevels
  'd pptn'
  'draw title \\\ GPCP precipitation "'date'" (mm d`a-1`n)'
  'run cbar 0.7'
  'set strsiz 0.1'
  'draw string 1.0 0.5 'min_max
  'set vpage  0 8.5 0 5.5'
  'set parea  1 7   1 4.5'
  'set grads off'
  'set gxout shaded'
  'set clevs 'clevels
  'd pptn5a'
  'draw title \\\ GPCP 5-month centered average precipitation (mm d`a-1`n) \ ("Method 1")'
  'run cbar 0.7'
  'set gxout stat'
  'd pptn5a'
  min_max = sublin(result,8)
  'set strsiz 0.1'
  'draw string 1.0 0.5 'min_max
  'print'
  pause()
  tt=tt+5
  n=n+1
endwhile

*------------------------------------------------
*  Method 2:  define only at desired time levels,
*  plot 5-month averages, difference with those
*  defined in "Method 1"...
*------------------------------------------------

tt = ts ; n = 1
while ( tt <= tf )
  say
  say 'Plot number 'n
  say '...time levels 'tt-2' 'tt+2
  'c'
  'set t 'tt
  'define pptn5b = ave(pptn,t-2,t+2)'
  'q time'
  date = subwrd(result,3) ; month = substr(date,6,3)
*
*  define contour levels
*
  'set gxout stat'
  'd pptn5b'
  min_max = sublin(result,8)
  rec14 = sublin(result,14)
  min = subwrd(min_max,4) ; max = subwrd(min_max,5)
  nmin = math_nint(min+0.5) ; nmax = math_nint(max-0.5)
*  rmin = min+0.5 ; rmax = max-0.5
  say ' min = 'min
  say ' max = 'max
  say ' nmin = 'nmin
  say ' nmax = 'nmax
*  say ' rmin = 'rmin
*  say ' rmax = 'rmax
  pmin = nmin
  pmax = nmax
  ncols = 10 ; nlevs = ncols+1
  int = (pmax-pmin)/ncols ; say 'int = 'int
  clevels = ''
  l=1
  clev = pmin
  while ( l <= nlevs )
    clevels = ''clevels' 'clev
    clev = clev + int
    l=l+1
  endwhile
  say 'clevels = 'clevels
  'set vpage  0 8.5 5.5 11'
  'set parea  1 7   1 4.5'
  'set grads off'
  'set gxout shaded'
  'set clevs 'clevels
  'd pptn5b'
  'draw title \\\ GPCP 5-month centered average precipitation "'date'" (mm d`a-1`n) \ ("Method 2")'
  'run cbar 0.7'
  'set strsiz 0.1'
  'draw string 1.0 0.5 'min_max
  'set vpage  0 8.5 0 5.5'
  'set parea  1 7   1 4.5'
  'set grads off'
  'set gxout shaded'
  'set clevs 'clevels
  'd pptn5b-pptn5a'
  'draw title \\\ "Method 2" - "Method 1" (mm d`a-1`n)'
  'run cbar 0.7'
  'set gxout stat'
  'd pptn5b-pptn5a'
  min_max = sublin(result,8)
  'set strsiz 0.1'
  'draw string 1.0 0.5 'min_max
  'print'
  pause()
  tt=tt+5
  n=n+1
endwhile

'disable print'
