function rad (args)
                                                                                                                   
expid    = subwrd(args,1)
output   = subwrd(args,2)
                                                                                                                   
* Define Seasons to Process
* -------------------------
seasons  = ''
       k = 3
while( k > 0 )
    season = subwrd(args,k)
if( season = '' )
    k = -1
else
    seasons = seasons % ' ' % season
k = k+1
endif
endwhile
'uppercase 'seasons
          seasons = result
                                                                                                                   
                                                                                                                   
* Initialize
* ----------
"reinit"
'set display color white'
'set clab off'
'c'
 

* Get Radiation Variables
* -----------------------
'getvar RADLW RADIATION'
        qname.1 = subwrd(result,1)
        qfile.1 = subwrd(result,2)
        scale.1 = subwrd(result,3)
        expdsc  = subwrd(result,4)
'getvar RADSW RADIATION'
        qname.2 = subwrd(result,1)
        qfile.2 = subwrd(result,2)
        scale.2 = subwrd(result,3)
'getvar O3 CHEMISTRY'
        qname.3 = subwrd(result,1)
        qfile.3 = subwrd(result,2)
        scale.3 = subwrd(result,3)

'getenv "GEOSUTIL"'
         geosutil = result
                                                                                                                  

* Experiment Datasets
* -------------------
'set dfile 'qfile.1
'setx'
'sety'
'setz'
'setdates'
'run getenv "BEGDATE"'
             begdate  = result
'run getenv "ENDDATE"'
             enddate  = result
'sett'
'makezf 'qname.1'*'scale.1' radlw z'
'seasonal radlwz'

'set dfile 'qfile.2
'setx'
'sety'
'setz'
'sett'
'makezf 'qname.2'*'scale.2' radsw z'
'seasonal radswz'

'set dfile 'qfile.3
'setx'
'sety'
'setz'
'sett'
'makezf 'qname.3'*'scale.3' o3 z'
'seasonal o3z'

'set dfile 'qfile.1
'getinfo zdim'
         zdim = result
'set z  'zdim
'getinfo level'
         ptop = result
'sety'
'setz'

* Loop over Seasons to Process
* ----------------------------
       m = 1
while( m > 0 )
    season = subwrd(seasons,m)
if( season = '' )
         m = -1
else
         m = m+1
         say 'Processing Season: 'season
                                                                                                                   
'count "'season'" 'begdate' 'enddate
  num = result

'run  'geosutil'/plots/rad/rad_1.gs 'expid' 'season' 'output' 'num' 'ptop
*'!    'geosutil'/plots/grads_util/make_gif -name radz  -rotate 90 -season 'season' -output 'output' &'
*'!    'geosutil'/plots/grads_util/make_gif -name radz2 -rotate 90 -season 'season' -output 'output' &'
*'!    'geosutil'/plots/grads_util/make_gif -name ozone -rotate 90 -season 'season' -output 'output' &'

endif
endwhile