* +-======-+ 
*  Copyright (c) 2003-2007 United States Government as represented by 
*  the Admistrator of the National Aeronautics and Space Administration.  
*  All Rights Reserved.
*  
*  THIS OPEN  SOURCE  AGREEMENT  ("AGREEMENT") DEFINES  THE  RIGHTS  OF USE,
*  REPRODUCTION,  DISTRIBUTION,  MODIFICATION AND REDISTRIBUTION OF CERTAIN 
*  COMPUTER SOFTWARE ORIGINALLY RELEASED BY THE UNITED STATES GOVERNMENT AS 
*  REPRESENTED BY THE GOVERNMENT AGENCY LISTED BELOW ("GOVERNMENT AGENCY").  
*  THE UNITED STATES GOVERNMENT, AS REPRESENTED BY GOVERNMENT AGENCY, IS AN 
*  INTENDED  THIRD-PARTY  BENEFICIARY  OF  ALL  SUBSEQUENT DISTRIBUTIONS OR 
*  REDISTRIBUTIONS  OF THE  SUBJECT  SOFTWARE.  ANYONE WHO USES, REPRODUCES, 
*  DISTRIBUTES, MODIFIES  OR REDISTRIBUTES THE SUBJECT SOFTWARE, AS DEFINED 
*  HEREIN, OR ANY PART THEREOF,  IS,  BY THAT ACTION, ACCEPTING IN FULL THE 
*  RESPONSIBILITIES AND OBLIGATIONS CONTAINED IN THIS AGREEMENT.
*  
*  Government Agency: National Aeronautics and Space Administration
*  Government Agency Original Software Designation: GSC-15354-1
*  Government Agency Original Software Title:  GEOS-5 GCM Modeling Software
*  User Registration Requested.  Please Visit http://opensource.gsfc.nasa.gov
*  Government Agency Point of Contact for Original Software:  
*  			Dale Hithon, SRA Assistant, (301) 286-2691
*  
* +-======-+ 
****************************************************************************
*                                                                          *
* This script computes the Seasonal Means and Anomalies of the INPUT field *
*                                                                          *
****************************************************************************
 
function anomaly (args)
name   = subwrd(args,1)
file   = subwrd(args,2)
offset = subwrd(args,3)


if( offset="" )
offset = 0
endif
if( file="_" )
file = ""
endif

'set t 1'
tmax = gettmax()
'getinfo time'
time = result
tbeg = gettbeg()

'getinfo xdim'
xdim = result
'getinfo ydim'
ydim = result
'set x 1 'xdim
'set y 1 'ydim

if( xdim>1 )
'set lon 0 360'
endif

tmam1 = 03-tbeg+1
tjja1 = 06-tbeg+1
tson1 = 09-tbeg+1
tdjf1 = 12-tbeg+1

if( tmam1<0 )
tmam1 = tmam1 + 12
endif

if( tjja1<0 )
tjja1 = tjja1 + 12
endif

if( tson1<0 )
tson1 = tson1 + 12
endif

if( tdjf1<0 )
tdjf1 = tdjf1 + 12
endif

tmam2 = tmam1 + 2
tjja2 = tjja1 + 2
tson2 = tson1 + 2
tdjf2 = tdjf1 + 2

say " "
say "Processing Seasonal Means for: "name
say "t_beg = "tbeg"  offset = "offset
say "t_max = "tmax
say " "
say "MAM: "tmam1" "tmam2
say "JJA: "tjja1" "tjja2
say "SON: "tson1" "tson2
say "DJF: "tdjf1" "tdjf2

* Define Climatology Field
* ------------------------
t1 = 1  + offset
t2 = 12 + offset

'set t 't1' 't2
'define 'name'clim'file' = ave('name',t+0,t='tmax',1yr)'
'modify 'name'clim'file' seasonal'
say " "
say "Defined "name"clim"file


* Define Anomolay Field
* ---------------------
'set t 1 'tmax
'define 'name'anom'file' = 'name' - 'name'clim'file
say "Defined "name"anom"file
'define 'name'var'file' = 'name'anom'file'*'name'anom'file
t1 = 1  + offset
t2 = 12 + offset

'set t 't1' 't2
'define 'name'varclim'file' = ave('name'var'file',t+0,t='tmax',1yr)'
'modify 'name'varclim'file' seasonal'
say " "
say "Defined "name"varclim"file


'set t 1'

'define 'name'djf'file' = ave('name'clim'file',t='tdjf1',t='tdjf2')'
say "Defined "name"djf"file

'define 'name'jja'file' = ave('name'clim'file',t='tjja1',t='tjja2')'
say "Defined "name"jja"file

'define 'name'son'file' = ave('name'clim'file',t='tson1',t='tson2')'
say "Defined "name"son"file

'define 'name'mam'file' = ave('name'clim'file',t='tmam1',t='tmam2')'
say "Defined "name"mam"file

'define 'name'ann'file' = ave('name'clim'file',t='t1',t='t2')'
say "Defined "name"ann"file
say " "

'define 'name'djfstd'file' = ave(sqrt('name'varclim'file'),t='tdjf1',t='tdjf2')'
say "Defined "name"djfstd"file

'define 'name'jjastd'file' = ave(sqrt('name'varclim'file'),t='tjja1',t='tjja2')'
say "Defined "name"jjastd"file

'define 'name'sonstd'file' = ave(sqrt('name'varclim'file'),t='tson1',t='tson2')'
say "Defined "name"sonstd"file

'define 'name'mamstd'file' = ave(sqrt('name'varclim'file'),t='tmam1',t='tmam2')'
say "Defined "name"mamstd"file

'define 'name'annstd'file' = ave(sqrt('name'varclim'file'),t='t1',t='t2')'
say "Defined "name"annstd"file

tjan = tdjf1 + 1
tfeb = tdjf1 + 2
tmar = tdjf1 + 3
tapr = tdjf1 + 4
tmay = tdjf1 + 5
tjun = tdjf1 + 6
tjul = tdjf1 + 7
taug = tdjf1 + 8
tsep = tdjf1 + 9
toct = tdjf1 + 10
tnov = tdjf1 + 11
tdec = tdjf1 + 12

'define 'name'jan'file' = 'name'clim'file'(t='tjan')'
'define 'name'feb'file' = 'name'clim'file'(t='tfeb')'
'define 'name'mar'file' = 'name'clim'file'(t='tmar')'
'define 'name'apr'file' = 'name'clim'file'(t='tapr')'
'define 'name'may'file' = 'name'clim'file'(t='tmay')'
'define 'name'jun'file' = 'name'clim'file'(t='tjun')'
'define 'name'jul'file' = 'name'clim'file'(t='tjul')'
'define 'name'aug'file' = 'name'clim'file'(t='taug')'
'define 'name'sep'file' = 'name'clim'file'(t='tsep')'
'define 'name'oct'file' = 'name'clim'file'(t='toct')'
'define 'name'nov'file' = 'name'clim'file'(t='tnov')'
'define 'name'dec'file' = 'name'clim'file'(t='tdec')'
say "Defined "name"jan-dec"file

'define 'name'janstd'file' = sqrt('name'varclim'file'(t='tjan'))'
'define 'name'febstd'file' = sqrt('name'varclim'file'(t='tfeb'))'
'define 'name'marstd'file' = sqrt('name'varclim'file'(t='tmar'))'
'define 'name'aprstd'file' = sqrt('name'varclim'file'(t='tapr'))'
'define 'name'maystd'file' = sqrt('name'varclim'file'(t='tmay'))'
'define 'name'junstd'file' = sqrt('name'varclim'file'(t='tjun'))'
'define 'name'julstd'file' = sqrt('name'varclim'file'(t='tjul'))'
'define 'name'augstd'file' = sqrt('name'varclim'file'(t='taug'))'
'define 'name'sepstd'file' = sqrt('name'varclim'file'(t='tsep'))'
'define 'name'octstd'file' = sqrt('name'varclim'file'(t='toct'))'
'define 'name'novstd'file' = sqrt('name'varclim'file'(t='tnov'))'
'define 'name'decstd'file' = sqrt('name'varclim'file'(t='tdec'))'
say "Defined "name"jan-dec std"file

  function gettmax(args)
  'query file '
  tmax  = sublin(result,5)
  tmax  = subwrd(tmax,12)
  return tmax

  function gettbeg(args)
  'getinfo month'
  month = result
  if(month="JAN")
  tbeg = 1
  endif
  if(month="FEB")
  tbeg = 2
  endif
  if(month="MAR")
  tbeg = 3
  endif
  if(month="APR")
  tbeg = 4
  endif
  if(month="MAY")
  tbeg = 5
  endif
  if(month="JUN")
  tbeg = 6
  endif
  if(month="JUL")
  tbeg = 7
  endif
  if(month="AUG")
  tbeg = 8
  endif
  if(month="SEP")
  tbeg = 9
  endif
  if(month="OCT")
  tbeg = 10
  endif
  if(month="NOV")
  tbeg = 11
  endif
  if(month="DEC")
  tbeg = 12
  endif
  return tbeg