* +-======-+ 
*  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 stats (args)
name = subwrd(args,1)
pi = 3.141592654

'getinfo xfreq'
         xfreq = result
'getinfo zfreq'
         zfreq = result
'getinfo file'
         curfile = result

'q dims'

xline = sublin(result,2)
lon1  = subwrd(xline,6)
lon2  = subwrd(xline,8)
yline = sublin(result,3)
lat1  = subwrd(yline,6)
lat2  = subwrd(yline,8)
zline = sublin(result,4)

if( zfreq = 'fixed' )
    zlev1 = subwrd(zline,9)
    zlev2 = zlev1
endif
if( zfreq = 'varying' )
    zlev1 = subwrd(zline,11)
    zlev2 = subwrd(zline,13)
endif

if( zfreq = 'fixed' & xfreq = 'varying' )
    'define avg = aave('name',       lon='lon1',lon='lon2',lat='lat1',lat='lat2')'
    'define rms = aave(pow('name',2),lon='lon1',lon='lon2',lat='lat1',lat='lat2')'
    'define var = ('name')*('name') - avg*avg'
    'define var = aave(var,          lon='lon1',lon='lon2',lat='lat1',lat='lat2')'
endif

if( xfreq = 'fixed' & zfreq = 'varying' )
         lon2  = lon1
         zmin  = zlev1
         zmax  = zlev2
     if( zlev2 < zlev1 )
         zmin  = zlev2
         zmax  = zlev1
      endif
      n = 0
      z = zmin
     'set z 'z
     'define avg = lat-lat'
     'define rms = lat-lat'
     while( z<=zmax )
        'set z 'z
        'define avg = avg + ave('name'        *cos(lat*'pi'/180),lat='lat1',lat='lat2')'
        'define rms = rms + ave( pow('name',2)*cos(lat*'pi'/180),lat='lat1',lat='lat2')'
         n = n + 1
         z = z + 1
     endwhile
        'define avg = avg / 'n
        'define rms = rms / 'n
        'set z 'zmin' 'zmax
        'define varz = ('name')*('name') - avg*avg'

      n = 0
      z = zmin
     'set z 'z
     'define var = lat-lat'
     while( z<=zmax )
        'set z 'z
        'define var = var + ave(varz*cos(lat*'pi'/180),lat='lat1',lat='lat2')'
         n = n + 1
         z = z + 1
     endwhile
        'define var = var / 'n
endif

'set lon 'lon1
'set lat 'lat1

'define std = sqrt(var)'
'define rms = sqrt(rms)'
'd avg'
   avg = subwrd(result,4)
'd std'
   std = subwrd(result,4)
'd rms'
   rms = subwrd(result,4)

'set lon 'lon1' 'lon2
'set lat 'lat1' 'lat2
'set z 'zlev1' 'zlev2

say 'Average of 'name' = 'avg' over lats: ('lat1','lat2')  lons: ('lon1','lon2')'
say 'Std.Dev of 'name' = 'std' over lats: ('lat1','lat2')  lons: ('lon1','lon2')'
say 'RtMnSqr of 'name' = 'rms' over lats: ('lat1','lat2')  lons: ('lon1','lon2')'

return avg' 'std' 'rms