*
*  Calculates the monthly precipitation climatology and anomalies, using
*  technique adapted from http://grads.iges.org/grads/gadoc/gradcomdmodify.html
*
*  To use, grads -p
*
'reinit'
'open /net/cjs/gpcp/pptn.ctl'
'enable print pptn.gx'
'set lon 0 360'
'set lat -90 90'
'set z 1'
*
*  define climatology
*
ts = 13 ; tsm1 = ts-1
tf = 120
'set t 1 12'
'define pptnclim = ave(pptn,t+'tsm1',t='tf',12)'
'modify pptnclim seasonal'
*
*  define anomaly
*
'set t 'ts' 'tf
'define pptnanom = pptn - pptnclim'
*
*  ---  results for April ---
*
tp = ts+3
'set t 'tp
'q time'
date = subwrd(result,3)
month = substr(date,6,3)
*
*  climatology
*
*  define contour levels
*
'set gxout stat'
'd pptnclim'
min_max = sublin(result,8)
rec14 = sublin(result,14)
min = subwrd(min_max,4) ; say 'min = 'min
max = subwrd(min_max,5) ; say 'max = 'max
nmin = math_nint(min+0.5) ; say 'nmin = 'nmin
nmax = math_nint(max-0.5) ; say 'nmax = 'nmax
rmin = min+0.5 ; say 'rmin = 'rmin
rmax = max-0.5 ; say 'rmax = 'rmax
pmin = nmin
pmax = nmax
ncols = 10 ; nlevs = ncols+1
int = (pmax-pmin)/ncols ; say 'int = 'int
clevels = ''
n=1
  clev = pmin
while ( n <= nlevs )
  clevels = ''clevels' 'clev
  clev = clev + int
  n=n+1
endwhile
say 'clevels = 'clevels
sigma_nm1 = subwrd(rec14,2)
var_nm1 = subwrd(rec14,3)
say 'sigma_nm1 = 'sigma_nm1
say 'var_nm1 = 'var_nm1

*
*  plot climatology
*
'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 pptnclim'
'draw title \\ GPCP Climatological Average Precipitation "'month'" (mm d`a-1`n)'
'run cbar 0.7'
'set strsiz 0.08'
'draw string 1.0 0.5 'min_max
*
*  plot anomaly
*
'set vpage  0 8.5 0 5.5'
'set parea  1 7   1 4.5'
'set grads off'
'set gxout stat'
'd pptnanom'
min_max = sublin(result,8)
'set gxout shaded'
'd pptnanom'
'draw title \\ GPCP Precipitation Anomaly 'date' (mm d`a-1`n)'
'run cbar 0.7'
'draw string 1.0 0.5 'min_max
'print'
'disable print'
