function epd (args)

cs=subwrd(args,1)
if (cs='-h')
  say ' '
  say ' Plots the Adimensional Available Potential Energy for the Horizontal Modes'
  say ' (Rossby + Mixed, Gravity West and Grevity East + Kelvin)'
  say ' in Function of the Zonal Wave Number'
  say ' '
  say ' Usage:'
  say ' epd cs var [Mmax]'
  say ' '
  say ' default: cs = br ; var = rb ; Mmax = 126'
  say ' '
  say ' Must be runned under grads -p (portrait)'
  say ' '
  say ' args meaning:'
  say ' '
  say ' cs : br (barotropic mode) or bc baroclinic mode closest to 250 m'
  say ' var : rb (Rossby + Mixed), gw (gravity west) or ge (gravity east + Kelvin)'
  say ' Mmax : Rhomb+1 (even), Rhomb is the zonal wave truncation (Rhomb odd number)'
  say '        is the total number of zonal waves plus the zonal mean (+1)'
  say ' '
  return
endif
if (cs='') ; cs='br' ; endif
if (cs!='br' & cs!='bc')
  say 'Wrong cs: 'cs' ; must be br or bc'
  return
endif
var=subwrd(args,2)
if (var='') ; var='rb' ; endif
if (var!='rb' & var!='gw' & var!='ge')
  say 'Wrong var: 'var' ; must be rb or gw or ge'
  return
endif
Mmax=subwrd(args,3)
if (Mmax='') ; Mmax=126 ; endif
mrb=Mmax-1
if (mrb < 100)
  Rhomb='RB00'mrb
endif
if (mrb >= 100 & mrb < 1000)
  Rhomb='RB0'mrb
endif
if (mrb >= 1000)
  Rhomb='RB'mrb
endif

rec=chkmdpt()
if (rec!=0)
  return
endif

* Getting Kmax (number of pressure levels) and csvm (reference case)
'run ../../Config/Case_Dates.gs 'MONAN' 1 'ZCAS' 'no
filecsdt='./Case_Dates.dgs'
rec=read(filecsdt)
lrec=sublin(rec,2)
Kmax=subwrd(lrec,1)
csvm=subwrd(lrec,5)
rec=close(filecsdt)

say ' '
say ' csvm = 'csvm'  Kmax = 'Kmax'  cs = 'cs'  var = 'var'  Mmax = 'Mmax
say ' '

if (cs='br')
  hh=01
  hl=00
endif
if (cs='bc')
  hh=08
  hl=07
endif
Trunc=Rhomb'L0'Kmax

if (var='rb') ; tit='Rossby and Mixed' ; endif
if (var='gw') ; tit='Gravity West' ; endif
if (var='ge') ; tit='Gravity East and Kelvin' ; endif

fcg='dir.cfg'
rec=read(fcg)
lrec=sublin(rec,2)
diri=subwrd(lrec,1)
rec=read(fcg)
lrec=sublin(rec,2)
dirfig=subwrd(lrec,1)%'png_'csvm'/'
rec=close(fcg)

filei='HF_'Trunc'_'csvm'_h0'hh'.ctl'

say ' '
say ' filei = 'filei
say ' '

'reinit'

dfilei=diri''filei
rec=chkopfl(dfilei)
if (rec!=0)
  return
endif

'open 'dfilei

if (cs='br')
  if (var='rb') ; clevs='0.0025 0.005 0.01 0.025 0.05 0.1 0.15 0.2 0.25' ; endif
  if (var='gw') ; clevs='0.12 0.15 0.20 0.235 0.245 0.2475 0.2485 0.249'; endif
  if (var='ge') ; clevs='0.175 0.20 0.235 0.245 0.248 0.249 0.25 0.26' ; endif
endif
if (cs='bc')
  if (var='rb') ; clevs='0.06 0.075 0.1 0.15 0.20 0.25 0.30' ; endif
  if (var='gw') ; clevs='0.1 0.125 0.15 0.175 0.20 0.21 0.22 0.23' ; endif
  if (var='ge') ; clevs='0.12 0.15 0.175 0.20 0.21 0.22 0.23 0.235 0.24 0.2425' ; endif
endif

'query ctlinfo'
k=1
while (k<=100)
  leh=sublin(result,k)
  if (subwrd(leh,1)='title')
    eps=subwrd(leh,5)
    heq=subwrd(leh,7)
    break
  endif
  k=k+1
endwhile
if (heq >= 1000) ; fmt='%-6.1f' ; endif
if (heq >= 100 & heq < 1000) ; fmt='%-6.2f' ; endif
if (heq >= 10 & heq < 100) ; fmt='%-6.3f' ; endif
if (heq <  10) ; fmt='%-6.4f' ; endif
heqi=math_format(fmt,heq)
say ' 'eps' 'heqi
figname='epd_'var'_'Trunc'_'csvm'_h0'hh

'set display grey white'
'clear'

'set grads off'
'set mproj off'

'set parea 1 8 1 9.5'
'set x 1 17'
'set y 1 24'
'set xlopts 1 8 0.175'
'set ylopts 1 8 0.175'
'set clopts 1 8 0.15'
'set ccolor 1'
'set cthick 8'
'set clab forced'
'set clevs 'clevs
'display epd'var
'draw title `1Potential Energy 'tit':\`2E'hl'`1 = 'eps'   `2H'hl'`1 = 'heq' m - 'csvm'\`1Contours: 'clevs
'draw ylab `1Meridional Index `2l'
'draw xlab `1Zonal Wave Number `2s'

'query gxinfo'
lx=sublin(result,3)
xa=subwrd(lx,4)
xb=subwrd(lx,6)
ly=sublin(result,4)
ya=subwrd(ly,4)
yb=subwrd(ly,6)

xs=xb
ys=ya-0.6
xrb=xs+0.1
yrb=ys+0.25
xra=xrb-1.8
yra=yrb-0.35
'set line 0 1 6'
'draw recf 'xra' 'yra' 'xrb' 'yrb
'set line 1 1 6'
'draw rec 'xra' 'yra' 'xrb' 'yrb

'set strsiz 0.16'
'set string 1 br 6 0'
'draw string 'xs' 'ys' `1'Trunc

fig=dirfig''figname'.png'
say ' '
say ' 'figname
say ' '
'printim 'fig' png white'

return

function chkmdpt()
'query gxinfo'
lm=sublin(result,2)
xb=subwrd(lm,4)
yb=subwrd(lm,6)
rec=0
if (yb < xb)
  say ' '
  say ' Error: must use grads -p (portrait)'
  say ' '
  rec=1
endif
return rec

function chkopfl(file)
res=read(file)
rec=sublin(res,1)
if (rec=0)
   lrec=sublin(res,2)
   dst=subwrd(lrec,1)
   if (dst!='dset') ; rec=1 ; endif
endif
if (rec!=0)
  say ' '
  say ' Error : 'file' does not exist or is not a ctl file'
  say ' '
endif
if (rec=0)
  res=close(file)
endif
return rec