function main(args)
*  public domain 5/2021                   Wesley Ebisuzaki
*
*  calculate 1991-2020 climatology
*  calculate anomaly
*
*  The data needs to be monthly means
*
* usage: anom.gs (anomaly) (grads expression)
*
*       example 1: anom.gs u ugrdprs
*
*       example 2: set lat 0
*                  set lon 0
*                  set time 1 last
*                  set lev 50
*                  anom uqbo aave(ugrdprs,lon=0,lon=360,lat=-10,lat=10(
*
* output:  (anomaly) is anomaly from 1991-2020 climatology
*          clim is 1991-2020 climatology
*
*         if anomaly is '' or "", only the climatolgy is computed
*
* requires: GrADS
*
anom=subwrd(args,1)
gexpr=subwrd(args,2)

if (anom = "''" | anom = '""') 
  say 'computes:  clim = 1991-2020 climatology'
else
  say 'computes: ' anom '=anomaly('gexpr')'
  say '           clim = 1991-2020 climatology'
endif

* Get the time environment
'q dim'
diminfo = result
line5 = sublin(diminfo,5)
if (subwrd(line5,3) = 'fixed')
  time1 = subwrd(line5,9)
  time2 = ''
else
  time1 = subwrd(line5,11)
  time2 = subwrd(line5,13)
endif

* calculate climo
'set time jan1991 dec1991'
'define clim=ave('gexpr',t+0,time=dec2020,+12)'
'modify clim seasonal'

* restore time environment
'set t ' time1 ' ' time2

* define anomaly
if (anom = "''" | anom = '""') 
  say 'only clim computed'
else
  say 'compute clim and ' anom
  'define 'anom'='gexpr'-clim'
endif

exit 0