* +-======-+ 
*  Copyright (c) 2003-2018 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 zplot (args)

expid   = subwrd(args,1)
export  = subwrd(args,2)
GC      = subwrd(args,3)
field   = subwrd(args,4)
season  = subwrd(args,5)
output  = subwrd(args,6)
logflag = subwrd(args,7)
nseason = subwrd(args,8)
bdate   = subwrd(args,9)
edate   = subwrd(args,10)
expdsc  = subwrd(args,11)
name    = subwrd(args,12)

say "export  = "export
say "LogFlag = "logflag

   blak = 1
undef   = 1e15
qmin    = undef
qmax    = undef
qcnt    = undef
cint    = undef
qtype   = "continuous"
prefix  = "0"
suffix  = " "

* Get Plotting Values from Resource File
* --------------------------------------
'run getenv "GEOSUTIL"'
             geosutil = result
PLOTRC = geosutil'/plots/grads_util/plot.rc'

say ''
                        'getresource 'PLOTRC' 'export'_'GC'_TITLE'  ; title   = result
                        'getresource 'PLOTRC' 'export'_'GC'_FACTOR' ; fact    = result

                        'getresource 'PLOTRC' 'export'_'GC'_Z_CCOLS'
if( result = 'NULL' ) ; 'getresource 'PLOTRC' 'export'_'GC'_CCOLS' ; endif
                                                            ccols = result

                        'getresource 'PLOTRC' 'export'_'GC'_Z_CLEVS'
if( result = 'NULL' ) ; 'getresource 'PLOTRC' 'export'_'GC'_CLEVS' ; endif
                                                            clevs = result

                        'getresource 'PLOTRC' 'export'_'GC'_CNTLEVS'; cntlevs = result
                        'getresource 'PLOTRC' 'export'_'GC'_YLAB'    ; ylab    = result
                        'getresource 'PLOTRC' 'export'_'GC'_YFMT'    ; yfmt    = result
                        'getresource 'PLOTRC' 'export'_'GC'_CBSCALE'; cbscale = result
                        'getresource 'PLOTRC' 'export'_'GC'_AXLIM'  ; axlim   = result
                        'getresource 'PLOTRC' 'export'_'GC'_ZLOG'   ; zlog    = result
say ''

if( zlog   != 'NULL' )
     logflag = zlog
endif

if( title   = 'NULL' )
  'getdesc 'name
     title  = result
endif
if( fact    = 'NULL' ) ; fact    = 1            ; endif
if( clevs   = 'NULL' ) ; clevs   = "0.0 0.0"    ; endif
if( ccols   = 'NULL' ) ; ccols   = "0.0 0.0"    ; endif
if( cbscale = 'NULL' ) ; cbscale =  0.65        ; endif
if( cntlevs = 'NULL' ) ; cntlevs = "0.0 0.0"    ; endif


*******************************************************************
**** Note:                                                     ****
****  prefix  = "0" will set lowest  shade to white            ****
****  suffix  = "0" will set highest shade to white            ****
****  contour = "0" will disable contours                      ****
****  clevs   = "lev1 lev2 ..." will use user-defined values   ****
****  qmin    = ...             will use user-defined value    ****
****  qmax    = ...             will use user-defined value    ****
****  qcnt    = ...             will use user-defined value    ****
****                                                           ****
****                                                           ****
****                                                           ****
****                                                           ****
****                                                           ****
****                                                           ****
*******************************************************************

'c'

* Temperature Tendencies
* ----------------------

if( export = 'TIM' )
     title = 'DTDT from Moist Physics (K/day)'
     fact  = 86400
    clevs  = '-5 -4 -3 -2 -1 -.8 -.6 -.4 -0.2 0.2 .4 .6 .8 1 2 3 4 5'
     qtype = 'plsmns'
endif

if( export = 'TIT' )
     title = 'DTDT from Turbulence (K/day)'
     fact  = 86400
    clevs  = '-1.8 -1.6 -1.4 -1.2 -1 -.8 -.6 -.4 -.2 1 2 3 4 5 6 7 8 9'
     qtype = 'plsmns'
endif


* Miscellaneous
* -------------
if( export = 'RI' & GC = 'TURBULENCE' )
     title = 'Richardson Number from Louis'
endif

if( export = 'SRFDIS' )
     title = 'Frictional Heating Rate from Surface Drag K/s Pa'
endif
if( export = 'RI'  & GC = 'MOIST' )
     title = 'Ice phase cloud particle effective radius  (x10`a-6`n m)'
     fact  = 1e6
endif
if( export = 'RL' )
     title = 'Liquid cloud particle effective radius  (x10`a-6`n m)'
     fact  = 1e6
endif
if( export = 'REV_CN' )
     title = 'Evaporation of Convective Rain  (g/kg/day)'
     fact  = 1000*86400
endif
if( export = 'REV_AN' )
     title = 'Evaporation of Anvil Rain  (g/kg/day)'
     fact  = 1000*86400
endif
if( export = 'REV_LS' )
     title = 'Evaporation of Large-Scale Rain  (g/kg/day)'
     fact  = 1000*86400
endif
if( export = 'RSU_AN' )
     title = 'Sublimation of Anvil Rain  (g/kg/day)'
     fact  = 1000*86400
endif
if( export = 'RSU_LS' )
     title = 'Sublimation of Large-Scale Rain  (g/kg/day)'
     fact  = 1000*86400
endif
if( export = 'RSU_CN' )
     title = 'Sublimation of Convective Rain  (g/kg/day)'
     fact  = 1000*86400
endif
if( export = 'TAUIR' )
     title = 'Longwave Optical Thickness at 800 cm`a-1`n'
     clevs = "5  10  15  20  25  30  35  40  45  50  60  70  80  90  100 120 120 130"
     ccols = "0 49 47 45 44 37 36 34 33 32 31 21 22 23 24 25 26 27 28"
endif
if( export = 'TAUCLI' )
     title = 'Shortwave Optical Thickness for Cloud Ice'
     clevs = "0.01  0.02  0.03  0.04  0.05  0.06   0.08  0.10  0.12  0.14  0.16  0.18  0.20  0.25  0.3 .35 .4 .45"
     ccols = "0 49 47 45 44 37 36 34 33 32 31 21 22 23 24 25 26 27 28"
endif
if( export = 'TAUCLW' )
     title = 'Shortwave Optical Thickness for Cloud Liquid'
     clevs = "0.1  0.2  0.3  0.4  0.5  0.6   0.8  1.0  1.2  1.4  1.6  1.8  2.0  2.5  3 3.5 4 4.5"
     ccols = "0 49 47 45 44 37 36 34 33 32 31 21 22 23 24 25 26 27 28"
endif


* Create Plot
*************

'getinfo zdim'
         zdim = result

'set display color white'
'set vpage off'
'set parea 1.5 10.5 0.75 7.5'

'set lat -90 90'
'set lon 0'
'define qmod = 'field''season'z'

'set grads off'
'set clab off'
if( blak = 1 ) ; 'black' ; endif
                                        'set zlog off'
if( logflag = "ON" | logflag = "on" | logflag = " ON" | logflag = " on" ) 
                                        'set zlog on'
                                        'setlevs' 
endif

if( axlim != "NULL" ) ; 'set axlim 'axlim ; endif
if( yfmt  != "NULL" ) ; 'set ylab  'yfmt  ; endif

'define  qdum = qmod*'fact

if( ccols = "0.0 0.0" )

  gg = xyz (xmin,xmax,ymin,ymax,zmin,zmax)
xmin = subwrd(gg,1)
xmax = subwrd(gg,2)
ymin = subwrd(gg,3)
ymax = subwrd(gg,4)
zmin = subwrd(gg,5)
zmax = subwrd(gg,6)

say ' XMIN: 'xmin'   XMAX: 'xmax
say ' YMIN: 'ymin'   YMAX: 'ymax
say ' ZMIN: 'zmin'   ZMAX: 'zmax

'set x 'xmin' 'xmax
'set y 'ymin' 'ymax

* Find QMIN and QMAX
* ------------------
     z  = zmin
'set z 'z

'minmax.simple qdum'
   qmax = subwrd(result,1)
   qmin = subwrd(result,2)
  qzmax = zmin
  qzmin = zmin
      z = z + 1

 while( z <= zmax )
'set z 'z
'minmax.simple qdum'
 dmax = subwrd(result,1)
 dmin = subwrd(result,2)
 if( dmax > qmax ) ; qmax = dmax ; qzmax = z ; endif
 if( dmin < qmin ) ; qmin = dmin ; qzmin = z ; endif
 z = z + 1
 endwhile
* ------------------

say ' QMIN: 'qmin'   QMAX: 'qmax
say 'QZMIN: 'qzmin'  QZMAX: 'qzmax

'd abs('qmin')'
        qmin = subwrd(result,4)
'd abs('qmax')'
        qmax = subwrd(result,4)

if( qmin > qmax ) 
    qmax = qmin
   'd abs('qzmin')'
           qzmin = subwrd(result,4)
   'd abs('qzmax')'
           qzmax = subwrd(result,4)
    if( qzmin > qzmax ) ; qzmax = qzmin ; endif
endif

* Determin if ZLOG is appropriate
* -------------------------------
if( zmin < zmax ) then
   'set lev 100'
   'getinfo zpos'
            zpos = result
   'set z 'zmin' 'zmax
   if( qzmax > zpos ) 
      'setlevs'
   else
      'set zlog off'
   endif
endif

'set z 'zmin' 'zmax
* -------------------------------

if( qmax > 0 )
   'd log10('qmax')'
    n = subwrd(result,4)
else
    n = 0
endif
say '    Log Factor: 'n

if( n<0 ) ; n = n-2 ; endif
'getint 'n
         n = result
if( n>0 )
    if( n<=2 ) 
        n = 0
     else
        n = n+2 
     endif
endif

if( n<0 )
    m = -n
else
    m =  n
endif

say 'Scaling Factor: 'n
if( n>0 )
    say 'define qdum = qdum / 1e'm 
        'define qdum = qdum / 1e'm 
else
    say 'define qdum = qdum * 1e'm 
        'define qdum = qdum * 1e'm 
endif

'minmax qdum'
   qmax = subwrd(result,1)
   qmin = subwrd(result,2)
say 'New MAX: 'qmax
say 'New MIN: 'qmin

     'shades  qdum  0'
      cint = result
else
    'set clevs 'clevs
    'set ccols 'ccols
         n = 0
endif

'set gxout shaded'
if( prefix = "0" | suffix = "0" ) ; 'black' ; endif
'd qdum'
'cbarn  -xmid 5.95 -ymid 0.25 -snum 'cbscale

'q shades'
num = sublin(result,1)
num = subwrd(num,5)
clevs = sublin(result,2)
clevs = subwrd(clevs,3)
k = 2
while( k<= num )
loc = k+1
line = sublin(result,loc)
clevs = clevs' 'subwrd(line,2)
k = k+1
endwhile

if( prefix = "0" | suffix = "0" ) ; 'black' ; endif
if( blak = 1 )                    ; 'black' ; endif
'set gxout contour'
'set clevs 'clevs
'set ccolor  1'
'd qdum'

if( zdim > 1 )       ; 'draw ylab Pressure (mb)' ; endif
if( ylab != "NULL" ) ; 'draw ylab 'ylab          ; endif

'set string 1 c 6 0'
'set strsiz .12'
'draw string 5.96 8.2 EXPID: 'expid'  'expdsc
if( n != 0 )
'draw string 5.96 7.8 'export':'GC'  'title'  'season' ('nseason')    (x10`a'n'`n)'
else
'draw string 5.96 7.8 'export':'GC'  'title'  'season' ('nseason')'
endif

* Print Beginning and Ending Dates
* --------------------------------
                date = getdate (bdate)
bmnthm = subwrd(date,1)
byearm = subwrd(date,2)
                date = getdate (edate)
emnthm = subwrd(date,1)
eyearm = subwrd(date,2)

'set string 1 l 4'
'set strsiz .08'
'draw string 0.50  8.20 Beg: 'bmnthm' 'byearm
'draw string 0.50  8.05 End: 'emnthm' 'eyearm
'set string 1 c 6'
* --------------------------------

'set string 1 c 6 0'

'myprint -name 'output'/hdiag_'export'.'GC'_z.'season

'set vpage off'
'set parea off'
return

function xyz (xmin,xmax,ymin,ymax,zmin,zmax)
'getinfo xfreq'
         xfreq = result
     if( xfreq = 'varying' )
         'getinfo xmin'
                  xmin = result
         'getinfo xmax'
                  xmax = result
     endif
     if( xfreq = 'fixed' )
         'getinfo xpos' 
                  xmin = result
                  xmax = result
     endif 

'getinfo yfreq'
         yfreq = result
     if( yfreq = 'varying' )
         'getinfo ymin'
                  ymin = result
         'getinfo ymax'
                  ymax = result
     endif
     if( yfreq = 'fixed' )
         'getinfo ypos' 
                  ymin = result
                  ymax = result
     endif 

'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 

return xmin' 'xmax' 'ymin' 'ymax' 'zmin' 'zmax

function getdate (date,month,year)
       num = 1
       bit = substr(date,num,1)
while( bit != '' )
       num = num+1
       bit = substr(date,num,1)
endwhile
       loc = num-7
     month = substr(date,loc  ,3)
      year = substr(date,loc+3,4)
return month' 'year