function zquad (args)

*******************************************************
****                 INPUT Variables               ****
*******************************************************

'numargs  'args
 numargs = result

DEBUG  = FALSE

        n   = 0
        num = 0
while ( num < numargs )
        num = num + 1

if( subwrd(args,num) = '-EXPID'  ) ; EXPID  = subwrd(args,num+1) ; endif
if( subwrd(args,num) = '-OUTPUT' ) ; OUTPUT = subwrd(args,num+1) ; endif
if( subwrd(args,num) = '-DEBUG'  ) ; DEBUG  = subwrd(args,num+1) ; endif

* Read SEASONS
* -----------
if( subwrd(args,num) = '-SEASON' )
     seasons = ''
           k = 1
    while( k > 0 )
           L = num + k
        season  = subwrd(args,L)
    if( season  = '' )
        k = -1
    else
        bit = substr(season,1,1)
        if( bit = '-' )
              k = -1
        else
              seasons = seasons % ' ' % season
              k = k+1
        endif
    endif
    endwhile
endif

endwhile
'uppercase 'seasons
            seasons = result
'numargs  ' seasons
 numseas =  result



* Initialize
* ----------
"reinit"
'set display color white'
'set csmooth on'
'c'


* Get Environment Variables
* -------------------------

'getenv "GEOSUTIL"'
         geosutil = result
'getenv "VERIFICATION"'
         verification = result


* Get Model Primary Fields
* ------------------------
'getvar PS DYN'
      pmname  = subwrd(result,1)
      pmfile  = subwrd(result,2)
      pmscale = subwrd(result,3)


'getvar U DYN'
say 'Result: 'result
      umname  = subwrd(result,1)
      umfile  = subwrd(result,2)
      umscale = subwrd(result,3)
       mdesc  = subwrd(result,4)

'getvar V DYN'
      vmname  = subwrd(result,1)
      vmfile  = subwrd(result,2)
      vmscale = subwrd(result,3)

'getvar T DYN'
      tmname  = subwrd(result,1)
      tmfile  = subwrd(result,2)
      tmscale = subwrd(result,3)

'getvar ZLE DYN'
      hmname  = subwrd(result,1)
      hmfile  = subwrd(result,2)
      hmscale = subwrd(result,3)

'getvar OMEGA DYN'
      wmname  = subwrd(result,1)
      wmfile  = subwrd(result,2)
      wmscale = subwrd(result,3)

'getvar Q MOIST'
      qmname  = subwrd(result,1)
      qmfile  = subwrd(result,2)
      qmscale = subwrd(result,3)

'define umscale = 'umscale
'define vmscale = 'vmscale
'define tmscale = 'tmscale
'define hmscale = 'hmscale
'define wmscale = 'wmscale
'define qmscale = 'qmscale

* Set Model CoVariance Names
* --------------------------
 result = setname (umname,vmname,umfile,m) 
          uvmname  = subwrd(result,1)
          uvmfile  = subwrd(result,2)
          uvmscale = umscale*vmscale

 result = setname (vmname,tmname,vmfile,m)
          vtmname  = subwrd(result,1)
          vtmfile  = subwrd(result,2)
          vtmscale = vmscale*tmscale

 result = setname (vmname,qmname,vmfile,m)
          vqmname  = subwrd(result,1)
          vqmfile  = subwrd(result,2)
          vqmscale = vmscale*qmscale

 result = setname (wmname,tmname,wmfile,m)
          wtmname  = subwrd(result,1)
          wtmfile  = subwrd(result,2)
          wtmscale = wmscale*tmscale

 result = setname (wmname,qmname,wmfile,m)
          wqmname  = subwrd(result,1)
          wqmfile  = subwrd(result,2)
          wqmscale = wmscale*qmscale

'define uvmscale = 'uvmscale
'define vtmscale = 'vtmscale
'define vqmscale = 'vqmscale
'define wtmscale = 'wtmscale
'define wqmscale = 'wqmscale

**************************************************
****            Compute Quadratics            ****
**************************************************

* Compute Zonal Means of PS Field
* -------------------------------
   'set dfile 'pmfile
   'setlons'
   'setdates'
   'set lat -90 90'
   'getdates'
   'set z 1'
   'run   setenv "ZDIM" 1'
   'makezf  'pmname' psmod z'
   'run unsetenv "ZDIM"'
   'set lon 0'
   'seasonal psmodz'
   'undefine psmodz'

* Compute Quadratics of Remaining Fields
* --------------------------------------
   'set dfile 'umfile
   'setlons'
   'setdates'
   'run getenv "BEGDATE"'
                begdate  = result
   'run getenv "ENDDATE"'
                enddate  = result
   'set lat -90 90'
   'getdates'
   'setz'

   'makequadz -U 'umname' -V 'vmname' -T 'tmname' -Q 'qmname' -H 'hmname' -W 'wmname' -UV 'uvmname' -VT 'vtmname' -VQ 'vqmname' -WT 'wtmname' -WQ 'wqmname' -TAG m'

   'set lon 0'
   'define   KEsm = 0.5*( Var'umname'sm + Var'vmname'sm )'
   'define   KEpm = 0.5*( Var'umname'pm + Var'vmname'pm )'
   'seasonal KEsm' ; 'undefine KEsm'
   'seasonal KEpm' ; 'undefine KEpm'

   'seasonal Var'umname'sm' ; 'undefine Var'umname'sm'
   'seasonal Var'vmname'sm' ; 'undefine Var'vmname'sm'
   'seasonal Var'tmname'sm' ; 'undefine Var'tmname'sm'
   'seasonal Var'qmname'sm' ; 'undefine Var'qmname'sm'
   'seasonal Var'hmname'sm' ; 'undefine Var'hmname'sm'
   'seasonal Var'wmname'sm' ; 'undefine Var'wmname'sm'

   'seasonal Var'umname'pm' ; 'undefine Var'umname'pm'
   'seasonal Var'vmname'pm' ; 'undefine Var'vmname'pm'
   'seasonal Var'tmname'pm' ; 'undefine Var'tmname'pm'
   'seasonal Var'qmname'pm' ; 'undefine Var'qmname'pm'
   'seasonal Var'hmname'pm' ; 'undefine Var'hmname'pm'
   'seasonal Var'wmname'pm' ; 'undefine Var'wmname'pm'

   'seasonal Cov'umname''vmname'sm' ; 'undefine Cov'umname''vmname'sm'
   'seasonal Cov'vmname''tmname'sm' ; 'undefine Cov'vmname''tmname'sm'
   'seasonal Cov'vmname''qmname'sm' ; 'undefine Cov'vmname''qmname'sm'
   'seasonal Cov'wmname''tmname'sm' ; 'undefine Cov'wmname''tmname'sm'
   'seasonal Cov'wmname''qmname'sm' ; 'undefine Cov'wmname''qmname'sm'

   'seasonal Cov'umname''vmname'pm' ; 'undefine Cov'umname''vmname'pm'
   'seasonal Cov'vmname''tmname'pm' ; 'undefine Cov'vmname''tmname'pm'
   'seasonal Cov'vmname''qmname'pm' ; 'undefine Cov'vmname''qmname'pm'
   'seasonal Cov'wmname''tmname'pm' ; 'undefine Cov'wmname''tmname'pm'
   'seasonal Cov'wmname''qmname'pm' ; 'undefine Cov'wmname''qmname'pm'


*********************************************************************
****   Loop over Possible Verification Datasets for Comparison   ****
*********************************************************************

' getnumrc 'geosutil'/plots/quad'
     rcinfo = result
     numrc  = subwrd( rcinfo,1 )
       num  = 1
       cnt  = 0

while( num <= numrc )
        loc = num + 1
     rcfile = subwrd( rcinfo,loc )


* Get OBS Primary Fields
* ----------------------
'getobs PS DYN 'rcfile
        poname = subwrd(result,1)
        pofile = subwrd(result,2)

'getobs U DYN 'rcfile
        uoname  = subwrd(result,1)
        uofile  = subwrd(result,2)
        uoscale = subwrd(result,3)
         odesc  = subwrd(result,4)
         obsid  = subwrd(result,5)

'getobs V DYN 'rcfile
        voname  = subwrd(result,1)
        vofile  = subwrd(result,2)
        voscale = subwrd(result,3)

'getobs T DYN 'rcfile
        toname  = subwrd(result,1)
        tofile  = subwrd(result,2)
        toscale = subwrd(result,3)

'getobs ZLE DYN 'rcfile
        honame  = subwrd(result,1)
        hofile  = subwrd(result,2)
        hoscale = subwrd(result,3)

'getobs OMEGA DYN 'rcfile
        woname  = subwrd(result,1)
        wofile  = subwrd(result,2)
        woscale = subwrd(result,3)

'getobs Q MOIST 'rcfile
        qoname  = subwrd(result,1)
        qofile  = subwrd(result,2)
        qoscale = subwrd(result,3)

'define uoscale = 'uoscale
'define voscale = 'voscale
'define toscale = 'toscale
'define hoscale = 'hoscale
'define woscale = 'woscale
'define qoscale = 'qoscale


* Set Obs CoVariance Names
* ------------------------
 result = setname (uoname,voname,uofile,o)
          uvoname  = subwrd(result,1)
          uvofile  = subwrd(result,2)
          uvoscale = uoscale*voscale

 result = setname (voname,toname,vofile,o)
          vtoname  = subwrd(result,1)
          vtofile  = subwrd(result,2)
          vtoscale = voscale*toscale

 result = setname (voname,qoname,vofile,o)
          vqoname  = subwrd(result,1)
          vqofile  = subwrd(result,2)
          vqoscale = voscale*qoscale

 result = setname (woname,toname,wofile,o)
          wtoname  = subwrd(result,1)
          wtofile  = subwrd(result,2)
          wtoscale = woscale*toscale

 result = setname (woname,qoname,wofile,o)
          wqoname  = subwrd(result,1)
          wqofile  = subwrd(result,2)
          wqoscale = woscale*qoscale

'define uvoscale = 'uvoscale
'define vtoscale = 'vtoscale
'define vqoscale = 'vqoscale
'define wtoscale = 'wtoscale
'define wqoscale = 'wqoscale


* Compute Zonal Means of PS Field
* -------------------------------
if( poname != 'NULL' )
   'set dfile 'pofile
   'setlons'
   'set lat -90 90'
   'getdates'
   'getinfo   tmin'
              tmin = result
   'getinfo   tmax'
              tmax = result
   'set t    'tmin
   'getinfo   date'
              begdateo = result
   'set t    'tmax
   'getinfo   date'
              enddateo = result
   'set time 'begdateo' 'enddateo
   'set z 1'
   'run   setenv "ZDIM" 1'
   'makezf  'poname' psobs z'
   'run unsetenv "ZDIM"'
   'set lon 0'
   'seasonal psobsz'
   'undefine psobsz'
else
   'set dfile 'uofile
   'setlons'
   'set lat -90 90'
   'getdates'
   'getinfo   tmin'
              tmin = result
   'getinfo   tmax'
              tmax = result
   'set t    'tmin
   'getinfo   date'
              begdateo = result
   'set t    'tmax
   'getinfo   date'
              enddateo = result
   'set time 'begdateo' 'enddateo
   'set z 1'
   'run   setenv "ZDIM" 1'
   'define psobs = 100000 + lon-lon + lat-lat'
   'makezf psobs psobs z'
   'run unsetenv "ZDIM"'
   'set lon 0'
   'seasonal psobsz'
   'undefine psobsz'
endif

   'run setenv   "BEGDATEO" 'begdateo
   'run setenv   "ENDDATEO" 'enddateo

* Compute Quadratics of Remaining Fields
* --------------------------------------
   'set dfile 'uofile
   'setlons'
   'set lat -90 90'
   'getdates'
   'getinfo   tmin'
              tmin = result
   'getinfo   tmax'
              tmax = result
   'set t    'tmin
   'getinfo   date'
              begdateo = result
   'set t    'tmax
   'getinfo   date'
              enddateo = result
   'set time 'begdateo' 'enddateo
   'setz'

   'makequadz -U 'uoname' -V 'voname' -T 'toname' -Q 'qoname' -H 'honame' -W 'woname' -UV 'uvoname' -VT 'vtoname' -VQ 'vqoname' -WT 'wtoname' -WQ 'wqoname' -TAG o'

   'set lon 0'
   'define   KEso = 0.5*( Var'uoname'so + Var'voname'so )'
   'define   KEpo = 0.5*( Var'uoname'po + Var'voname'po )'
   'seasonal KEso' ; 'undefine KEso'
   'seasonal KEpo' ; 'undefine KEpo'

   'seasonal Var'uoname'so' ; 'undefine Var'uoname'so'
   'seasonal Var'voname'so' ; 'undefine Var'voname'so'
   'seasonal Var'toname'so' ; 'undefine Var'toname'so'
   'seasonal Var'qoname'so' ; 'undefine Var'qoname'so'
   'seasonal Var'honame'so' ; 'undefine Var'honame'so'
   'seasonal Var'woname'so' ; 'undefine Var'woname'so'

   'seasonal Var'uoname'po' ; 'undefine Var'uoname'po'
   'seasonal Var'voname'po' ; 'undefine Var'voname'po'
   'seasonal Var'toname'po' ; 'undefine Var'toname'po'
   'seasonal Var'qoname'po' ; 'undefine Var'qoname'po'
   'seasonal Var'honame'po' ; 'undefine Var'honame'po'
   'seasonal Var'woname'po' ; 'undefine Var'woname'po'

   'seasonal Cov'uoname''voname'so' ; 'undefine Cov'uoname''voname'so'
   'seasonal Cov'voname''toname'so' ; 'undefine Cov'voname''toname'so'
   'seasonal Cov'voname''qoname'so' ; 'undefine Cov'voname''qoname'so'
   'seasonal Cov'woname''toname'so' ; 'undefine Cov'woname''toname'so'
   'seasonal Cov'woname''qoname'so' ; 'undefine Cov'woname''qoname'so'

   'seasonal Cov'uoname''voname'po' ; 'undefine Cov'uoname''voname'po'
   'seasonal Cov'voname''toname'po' ; 'undefine Cov'voname''toname'po'
   'seasonal Cov'voname''qoname'po' ; 'undefine Cov'voname''qoname'po'
   'seasonal Cov'woname''toname'po' ; 'undefine Cov'woname''toname'po'
   'seasonal Cov'woname''qoname'po' ; 'undefine Cov'woname''qoname'po'


*******************************************************
****         Loop over Seasons to Process          ****
*******************************************************

ptop.1 = 100
ptop.2 = 1

       m = 1
while( m > 0 )
    season = subwrd(seasons,m)
if( season = '' )
         m = -1
else
         m = m+1
         say 'Processing Season: 'season

'set dfile 'umfile
'count "'season'" 'begdate' 'enddate
 nmod =  result

'set dfile 'uofile
'count "'season'"'
 nobs =  result

* Make Variance Plots
* -------------------
          VARSM = umname' 'vmname' 'tmname' 'qmname' 'wmname' 'hmname
          VARSO = uoname' 'voname' 'toname' 'qoname' 'woname' 'honame
'numargs 'VARSM
          nvars = result

n = 1
while( n<=nvars )
       varmname = subwrd( VARSM,n )
       varoname = subwrd( VARSO,n )

if( varmname != 'NULL' & varoname != 'NULL' )

'uppercase 'varmname
            varname = 'Var_'result

'set dfile 'umfile
'setz'
'define qsqsmz = sqrt(Var'varmname'sm'season')'
'define qpqpmz = sqrt(Var'varmname'pm'season')'
'set dfile 'uofile
'setz'
'define qsqsoz = sqrt(Var'varoname'so'season')'
'define qpqpoz = sqrt(Var'varoname'po'season')'

k = 1
while( k<=2 )
ptop = ptop.k

           flag = ""
   while ( flag = "" )
  'run 'geosutil'/plots/quad/zplot -EXPID 'EXPID' -EXPORT 'varname' -SEASON 'season' -NMOD 'nmod' -NOBS 'nobs' -OUTPUT 'OUTPUT' -PFILE 'umfile' -OFILE 'uofile' -OBDATE 'begdateo' -OEDATE 'enddateo' -ODESC 'odesc' -ONAME 'obsid' -MDESC 'mdesc' -PTOP 'ptop
        if( DEBUG = "debug" )
            say "Hit  ENTER  to repeat plot or non-blank to continue"
            pull flag
        else
            flag = "next"
        endif
       'c'
   endwhile

k = k+1
endwhile

           flag = ""
   while ( flag = "" )
  'run 'geosutil'/plots/quad/lplot -EXPID 'EXPID' -EXPORT 'varname' -SEASON 'season' -NMOD 'nmod' -NOBS 'nobs' -OUTPUT 'OUTPUT' -PFILE 'umfile' -OFILE 'uofile' -OBDATE 'begdateo' -OEDATE 'enddateo' -ODESC 'odesc' -ONAME 'obsid' -PTOP 'ptop
        if( DEBUG = "debug" )
            say "Hit  ENTER  to repeat plot or non-blank to continue"
            pull flag
        else
            flag = "next"
        endif
       'c'
   endwhile

endif

n = n+1
endwhile

* Make KE Plots
* -------------
'set dfile 'umfile
'setz'
'define qsqsmz = KEsm'season
'define qpqpmz = KEpm'season
'set dfile 'uofile
'setz'
'define qsqsoz = KEso'season
'define qpqpoz = KEpo'season

k = 1
while( k<=2 )
ptop = ptop.k

           flag = ""
   while ( flag = "" )
  'run 'geosutil'/plots/quad/zplot -EXPID 'EXPID' -EXPORT Var_KE -SEASON 'season' -NMOD 'nmod' -NOBS 'nobs' -OUTPUT 'OUTPUT' -PFILE 'umfile' -OFILE 'uofile' -OBDATE 'begdateo' -OEDATE 'enddateo' -ODESC 'odesc' -ONAME 'obsid' -MDESC 'mdesc' -PTOP 'ptop
        if( DEBUG = "debug" )
            say "Hit  ENTER  to repeat plot or non-blank to continue"
            pull flag
        else
            flag = "next"
        endif
       'c'
   endwhile

k = k+1
endwhile

           flag = ""
   while ( flag = "" )
  'run 'geosutil'/plots/quad/lplot -EXPID 'EXPID' -EXPORT Var_KE -SEASON 'season' -NMOD 'nmod' -NOBS 'nobs' -OUTPUT 'OUTPUT' -PFILE 'umfile' -OFILE 'uofile' -OBDATE 'begdateo' -OEDATE 'enddateo' -ODESC 'odesc' -ONAME 'obsid' -PTOP 'ptop
        if( DEBUG = "debug" )
            say "Hit  ENTER  to repeat plot or non-blank to continue"
            pull flag
        else
            flag = "next"
        endif
       'c'
   endwhile


* Make CoVariance Plots
* ---------------------
                  nvars = 1
          COVARSM.nvars = umname' 'vmname ; nvars = nvars + 1
          COVARSM.nvars = vmname' 'tmname ; nvars = nvars + 1
          COVARSM.nvars = vmname' 'qmname ; nvars = nvars + 1
          COVARSM.nvars = wmname' 'tmname ; nvars = nvars + 1
          COVARSM.nvars = wmname' 'qmname ; nvars = nvars + 1
                  nvars = nvars - 1

                  nvars = 1
          COVARSO.nvars = uoname' 'voname ; nvars = nvars + 1
          COVARSO.nvars = voname' 'toname ; nvars = nvars + 1
          COVARSO.nvars = voname' 'qoname ; nvars = nvars + 1
          COVARSO.nvars = woname' 'toname ; nvars = nvars + 1
          COVARSO.nvars = woname' 'qoname ; nvars = nvars + 1
                  nvars = nvars - 1

n = 1
while( n<=nvars )
       varmname1 = subwrd( COVARSM.n,1 )
       varmname2 = subwrd( COVARSM.n,2 )
       varoname1 = subwrd( COVARSO.n,1 )
       varoname2 = subwrd( COVARSO.n,2 )

'uppercase 'varmname1
             varname1 = result
'uppercase 'varmname2
             varname2 = result
               cname  = 'Cov_'varname1'_'varname2

'set dfile 'umfile
'setz'
'define qsqsmz = Cov'varmname1''varmname2'sm'season
'define qpqpmz = Cov'varmname1''varmname2'pm'season
'set dfile 'uofile
'setz'
'define qsqsoz = Cov'varoname1''varoname2'so'season
'define qpqpoz = Cov'varoname1''varoname2'po'season

k = 1
while( k<=2 )
ptop = ptop.k

           flag = ""
   while ( flag = "" )
  'run 'geosutil'/plots/quad/zplot -EXPID 'EXPID' -EXPORT 'cname' -SEASON 'season' -NMOD 'nmod' -NOBS 'nobs' -OUTPUT 'OUTPUT' -PFILE 'umfile' -OFILE 'uofile' -OBDATE 'begdateo' -OEDATE 'enddateo' -ODESC 'odesc' -ONAME 'obsid' -MDESC 'mdesc' -PTOP 'ptop
        if( DEBUG = "debug" )
            say "Hit  ENTER  to repeat plot or non-blank to continue"
            pull flag
        else
            flag = "next"
        endif
       'c'
   endwhile

k = k+1
endwhile

           flag = ""
   while ( flag = "" )
  'run 'geosutil'/plots/quad/lplot -EXPID 'EXPID' -EXPORT 'cname' -SEASON 'season' -NMOD 'nmod' -NOBS 'nobs' -OUTPUT 'OUTPUT' -PFILE 'umfile' -OFILE 'uofile' -OBDATE 'begdateo' -OEDATE 'enddateo' -ODESC 'odesc' -ONAME 'obsid' -PTOP 'ptop
        if( DEBUG = "debug" )
            say "Hit  ENTER  to repeat plot or non-blank to continue"
            pull flag
        else 
            flag = "next"
        endif
       'c'
   endwhile

n = n+1
endwhile

endif
endwhile

num = num + 1
endwhile

return

function setname (name1,name2,cfile,tag)

    flag = false
    if( name1 != 'NULL' & name2 != 'NULL' )
        cname  = 'Cov_'name1'_'name2
       'lowercase 'cname
                  lcname = result
       'set dfile 'cfile
       'query file'
        numvar = sublin(result,6)
        numvar = subwrd(numvar,5)
        m = 1
        while ( m<numvar+1 )
        field = sublin(result,6+m)
        field = subwrd(field,1)
              if( lcname=field )
              flag = true
              endif
        m = m+1
        endwhile
    endif
    if (flag=false)
        cname = 'Cov'name1''name2''tag
    endif
return cname' 'cfile