* +-======-+ 
*  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
*  
* +-======-+ 
function statmak (args)
field  = subwrd  (args,1)
tag    = subwrd  (args,2)

'numargs  'args
 numargs = result

* Query Level Environment
* -----------------------
'getinfo zfreq'
         zfreq = result
     if( zfreq = 'varying' )
         'getinfo zmin'
                  zmin = result
         'getinfo zmax'
                  zmax = result
     endif
     if( zfreq = 'fixed' )
         'getinfo zpos'
          zmin = result
          zmax = result
     endif

'run getenv "GEOSUTIL"'
             geosutil = result

'getinfo numfiles'
         numfiles = result

'open 'geosutil'/plots/grads_util/lwmask1440721.tabl'
'getinfo numfiles'
         newfile = result

'set dfile 'newfile
'set t 1'
'set z 1'
'define mask = lwmask'
'close 'newfile

'set dfile 1'
'set t 1'
'setlons'
'sety'
'set z 'zmin' 'zmax


* Define Number of Forecast Days and Time Interval (hrs)
* ------------------------------------------------------
'run getenv "NDAY"'
             nday = result
'run getinfo tinc'
             tinc = result
if( nday = "NULL" ) ; nday = 5 ; endif
if( tinc = "NULL" ) ; tinc = 6 ; endif
                      nmax = 1 + nday*(24/tinc)

'getinfo tdim'
         tdim = result
         tbeg = 1-(nmax-tdim)
'set t ' tbeg ' 'tdim

* Compute forecast statistics
* ---------------------------
*  fma: forecast minus analysis
*  fmc: forecast minus climatology
*  mes: mean error squared
*  mse: mean square error
*  rms: root mean square error
*  std: standard deviation
* --------------------------------------------------
'define 'field'fmean = lat-lat+lon-lon'
'define 'field'amean = lat-lat+lon-lon'
'define 'field'cmean = lat-lat+lon-lon'
'define 'field'fma   = lat-lat+lon-lon'
'define 'field'fmc   = lat-lat+lon-lon'
'define 'field'mes   = lat-lat+lon-lon'
'define 'field'mse   = lat-lat+lon-lon'
'define 'field'rms   = lat-lat+lon-lon'
'define 'field'std   = lat-lat+lon-lon'
'define 'field'var   = lat-lat+lon-lon'

n = 1
while ( n <= numfiles )
say 'Processing Field: 'field' for File: 'n
'define 'field'fmean = 'field'fmean +     'field'f.'n
'define 'field'amean = 'field'amean +     'field'a.'n
'define 'field'cmean = 'field'cmean +     'field'c.'n
'define 'field'fma   = 'field'fma   +     'field'f.'n'-'field'a.'n
'define 'field'fmc   = 'field'fmc   +     'field'f.'n'-'field'c.'n
'define 'field'mse   = 'field'mse   + pow('field'f.'n'-'field'a.'n',2)'
n = n + 1
endwhile

'define 'field'fmean'tag' = 'field'fmean/'numfiles
'define 'field'amean'tag' = 'field'amean/'numfiles
'define 'field'cmean'tag' = 'field'cmean/'numfiles
'define 'field'fma'tag'   = 'field'fma  /'numfiles
'define 'field'fmc'tag'   = 'field'fmc  /'numfiles
'define 'field'mse'tag'   = 'field'mse  /'numfiles

n = 1
while ( n <= numfiles )
'define 'field'mes = 'field'mes + pow('field'fmean'tag'-'field'amean'tag',2)'
n = n + 1
endwhile

'define 'field'mes'tag' =      'field'mes/'numfiles
'define 'field'rms'tag' = sqrt('field'mse'tag')'
'define 'field'var'tag' =      'field'mse'tag'-'field'mes'tag
'define 'field'std'tag' = sqrt('field'mse'tag'-'field'mes'tag')'

if( zfreq = 'varying' )
'makez 'field'fmean'tag' z'
'makez 'field'amean'tag' z'
'makez 'field'cmean'tag' z'
'makez 'field'fma'tag'   z'
'makez 'field'fmc'tag'   z'
'makez 'field'mes'tag'   z'
'makez 'field'mse'tag'   z'
'define 'field'rms'tag'z = sqrt('field'mse'tag'z)'
'define 'field'var'tag'z =      'field'mse'tag'z-'field'mes'tag'z'
'define 'field'std'tag'z = sqrt('field'mse'tag'z-'field'mes'tag'z)'
endif

return