function xsec (args)

caso=subwrd(args,1)
if (caso='-h')
  say ' '
  say ' Plots a Vertical Cross-Section of the Vertical Mode Energy'
  say ' or Percentage for the Classes 0 to 4 individually'
  say ' '
  say ' Usage:'
  say ' xsec [caso] [epoca] [casestudy] [mean] [Mmax] [etc] [pt] [dp]'
  say ' '
  say ' Default: caso = MONAN ; epoca=1 ; casestudy=ZCAS ; mean=no ; Mmax=126 ; etc=etc1 ; pt=0 ; dp=Perc'
  say ' '
  say ' Must be runned under grads -p (portrait)'
  say ' '
  say ' args meaning:'
  say ' '
  say ' caso : five letters idenfying the case'
  say ' epoca : number idenfying the datei and datef'
  say ' casestudy : name of the case study'
  say ' mean : yes if the input is a time mean or no otherwise'
  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 ' etc : is the modal energy of the class n'
  say ' pt : top pressure of vertical cross-section'
  say '      If pt=0 as entry, pt will be set as the upper level of the csvm'
  say ' dp : is to display percentage if Perc or the modal energy if != Perc'
  say ' '
  return
endif
* caso : five letters input case identifier
if (caso='') ; caso='MONAN' ; endif
epoca=subwrd(args,2)
* epoca : number idenfying the datei and datef
if (epoca='') ; epoca=1 ; endif
* casestudy : name of the case study
casestudy=subwrd(args,3)
if (casestudy='') ; casestudy='ZCAS' ; endif
* mean : yes if the input is a time mean or no otherwise
mean=subwrd(args,4)
if (mean='') ; mean='no' ; endif
* Mmax is the Rhomboibal truncation plus 1
Mmax=subwrd(args,5)
if (Mmax='') ; Mmax=126 ; endif
etc=subwrd(args,6)
* etc is the modal energy of the class n
if (etc='') ; etc='etc1' ; endif
* If pt=0 as entry, pt will be set as the upper level of the case
pt=subwrd(args,7)
if (pt='') ; pt=0 ; endif
* dp is to display percentage if dp=Perc or the modal energy if dp!=Perc
dp=subwrd(args,8)
if (dp='') ; dp='Perc' ; endif
if (dp!='Perc') ; dp='Ener' ; endif
pn=substr(etc,4,1)
vn=substr(etc,1,3)

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

say ' '
say ' caso = 'caso'  epoca = 'epoca'  casestudy = 'casestudy'  mean = 'mean
say ' Mmax = 'Mmax'  etc = 'etc'  pt = 'pt'  dp = 'dp
say ' '

* Getting Kmax (number of pressure levels), datei (initial datei - icn),
* datef (forecast date - fct) (if initial condition or reanalysis datei=datef)
* and prev (preffix: P.icn. for initial condition and P.fct. for forecast)
* for the caso defined by epoca
'run ../../Config/Case_Dates.gs 'caso' 'epoca' 'casestudy' 'mean
filecsdt='./Case_Dates.dgs'
rec=read(filecsdt)
lrec=sublin(rec,2)
Kmax=subwrd(lrec,1)
datei=subwrd(lrec,2)
datef=subwrd(lrec,3)
prev=subwrd(lrec,4)
csvm=subwrd(lrec,5)
rec=read(filecsdt)
lrec=sublin(rec,2)
z0a=subwrd(lrec,1) ; z0b=subwrd(lrec,2)
z1a=subwrd(lrec,3) ; z1b=subwrd(lrec,4)
z2a=subwrd(lrec,5) ; z2b=subwrd(lrec,6)
z3a=subwrd(lrec,7) ; z3b=subwrd(lrec,8)
z4a=subwrd(lrec,9) ; z4b=subwrd(lrec,10)
rec=close(filecsdt)
'!rm -f 'filecsdt
say ' caso = 'caso'   Kmax = 'Kmax'    datei = 'datei'    datef = 'datef'   prev = 'prev'   csvm = 'csvm
say ' '
say ' z0a - 'z0a'   z0b = 'z0b
say ' z1a - 'z1a'   z1b = 'z1b
say ' z2a - 'z2a'   z2b = 'z2b
say ' z3a - 'z3a'   z3b = 'z3b
say ' z4a - 'z4a'   z4b = 'z4b

if (pt=0)
* Read Pressure Levels
  if (Kmax < 100)
    FilePres='../../Pressure_Levels/PressureLevels_'csvm'.L0'Kmax
  else
    FilePres='../../Pressure_Levels/PressureLevels_'csvm'.L'Kmax
  endif
  nRec=Kmax/10
  if (math_int(nRec) != nRec) ; nRec=math_int(nRec)+1 ; endif
  PrsLev=''
  n=1
  while (n <= nRec)
    rec=read(FilePres)
    lrec=sublin(rec,2)
    PrsLev=PrsLev''%''lrec
    n=n+1
  endwhile
  rec=close(FilePres)
  pt=subwrd(PrsLev,Kmax)
endif
say ' '
say ' pt = 'pt

* Getting output resolution: Rhomb
'run ../../Config/Rhomb_Resol.gs 'Mmax
fileres='./Rhomb_Resol.dgs'
rec=read(fileres)
lrec=sublin(rec,2)
rec=close(fileres)
'!rm -fr 'fileres
Rhomb=subwrd(lrec,1)
say ' '
say ' Mmax = 'Mmax'   Rhomb = 'Rhomb
if (Kmax < 100)
  _Trunc=Rhomb'L0'Kmax
else
  _Trunc=Rhomb'L'Kmax
endif

mh='jan feb mar apr may jun jul aug sep oct nov dec'
yyi=substr(datei,1,4)
mmi=substr(datei,5,2)
ddi=substr(datei,7,2)
hhi=substr(datei,9,2)
msi=subwrd(mh,mmi)
dategi=hhi'z'ddi''msi''yyi
yyf=substr(datef,1,4)
mmf=substr(datef,5,2)
ddf=substr(datef,7,2)
hhf=substr(datef,9,2)
msf=subwrd(mh,mmf)
dategf=hhf'z'ddf''msf''yyf

if (datei=datef)
  tit='`1Total Modal Energy Percantage Case: 'caso'-'_Trunc'\`1Analysis for Date 'dategf
  if (mean='yes')
    tit='`1Total Modal Energy Percantage Case: 'caso'-'_Trunc'\`1Analysis: Mean for 00Z23JAN2024 to 00Z29JAN2024'
  endif
else
  tit='`1Total Modal Energy Percantage Case: 'caso'-'_Trunc'\`1Forecast for Dates 'dategi' - 'dategf
  if (mean='yes')
    tit='`1Total Modal Energy Percantage Case: 'caso'-'_Trunc'\`13 Day Forecast: Mean for 00Z23JAN2024 to 00Z29JAN2024'
  endif
endif

datego=dategi'_'dategf
dateo=datei''datef
fcg='dir.cfg'
rec=read(fcg)
lrec=sublin(rec,1)
lrec=sublin(rec,2)
diri=subwrd(lrec,1)
rec=read(fcg)
lrec=sublin(rec,2)
dirh=subwrd(lrec,1)
rec=read(fcg)
lrec=sublin(rec,2)
_dirfig=subwrd(lrec,1)%'png_'caso'/'
rec=close(fcg)

filei='ENCL'caso''dateo''prev''_Trunc'.ctl'
fileh='Vertical_Functions_'csvm'.L0'Kmax'.ctl'

say ' '
say ' diri = 'diri
say ' filei = 'filei
say ' dirh = 'dirh
say ' fileh = 'fileh

* Getting longitudes and latitudes limites of the Cross-Section
'run ../../Config/Cross_Section.gs 'casestudy
filecrssct='./Cross_Section.dgs'
rec=read(filecrssct)
lrec=sublin(rec,2)
LonA=subwrd(lrec,1)
LonB=subwrd(lrec,2)
LatA=subwrd(lrec,3)
LatB=subwrd(lrec,4)
rec=close(filecrssct)
'!rm -fr 'filecrssct
say ' '
say ' LonA = 'LonA'    LonB = 'LonB
say ' LatA = 'LatA'    LatB = 'LatB
say ' '
if (LonA > 360.0)
  say ' Longitudes and Latitudes not defined'
  say ' for the Case Study : 'casestudy
  return
endif

'reinit'

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

'open 'dfilei
'open 'dfileh

'set display color white'
'clear'
'set grads off'
'set parea 1 8 1 9.5'

'set lon 'LonA' 'LonB
if (LatA < LatB)
'set lat 'LatA' 'LatB
else
'set lat 'LatB' 'LatA
endif
'set lev 1000 'pt
if (dp='Perc')
  'define var=100*'etc'/(etc0+etc1+etc2+etc3+etc4)'
else
  if (etc='ett')
    'define var=etc0+etc1+etc2+etc3+etc4'
  else
    'define var='etc
  endif
endif
'set zlog on '
'set ylevs 1000 850 750 600 500 400 300 200 150 100 75 50 30 20 10 7 5 3 2 1'
'set x 1'
'set y 1'
'collect 1 free'
if (LonA!=LonB)
dlon=(LonB-LonA)/25.0
lon=LonA
while (lon<=LonB)
  lat=LatA+(LatB-LatA)*(lon-LonA)/(LonB-LonA)
  'collect 1 gr2stn(var,'lon','lat')'
  lon=lon+dlon
endwhile
endif
if (LatA!=LatB)
dlat=(LatB-LatA)/25.0
lat=LatA
while (lat<=LatB)
  lon=LonA+(LonB-LonA)*(lat-LatA)/(LatB-LatA)
  'collect 1 gr2stn(var,'lon','lat')'
  lat=lat+dlat
endwhile
endif
'set lon 'LonA' 'LonB
'set warn off'
'set xaxis 'LonA' 'LonB
'set gxout contour' 
'set xlopts 0'
'set ylopts 1 8 0.15'
'set clopts 1 8 0.12'
'set clab forced'
'set cthick 8'
'display coll2gr(1,-u)'
pa=subwrd(result,2)
pb=subwrd(result,4)
ip=subwrd(result,6)
ct='Contours: 'pa' to 'pb' by 'ip

nd=10
'query gxinfo'
lnx=sublin(result,3)
lny=sublin(result,4)
x0=subwrd(lnx,4)
xm=subwrd(lnx,6)
y0=subwrd(lny,4)
ym=subwrd(lny,6)
dx=(xm-x0)/nd
dsx=(LonB-LonA)/nd
dsy=(LatB-LatA)/nd

xa=0.875
ya=0.8*y0
xb=xa
yb=0.55*y0
xc=0.5*(x0+xm)
yc=0.3*y0
xd=0.3*x0
yd=0.5*(y0+ym)
xlo=xa-0.4
xla=xlo
'set line 1 1 8'
'set strsiz 0.15'
'set string 1 c 8 0'
'draw string 'xlo' 'ya' `1Lon:'
'draw string 'xla' 'yb' `1Lat:'
'draw string 'xc' 'yc' `1Cross Section A to B : 'ct
'set string 1 c 8 90'
'draw string 'xd' 'yd' `1Pressure (hPa)'
'set string 1 c 8 0'

x=x0
y=ya
s=LonA
t=1
while (t<nd+2)
if (s>'180')
s=s-360
endif
s=math_nint(100*s)/100
if (s<'0')
sa=-s
'draw string 'x' 'y' `1'sa'W'
if (t>1 & t<nd+1)
   'draw line 'x' 'y0' 'x' 'y0+0.1
endif
if (t=1)
x1=sa'W'
endif
if (t=nd+1)
x2=sa'W'
endif
else
if (s!='0')
if (s='180')
'draw string 'x' 'y' `1'180
else
'draw string 'x' 'y' `1's'E'
endif
else
'draw string 'x' 'y' `1'0
endif
if (t=1)
x1=s'E'
endif
if (t=nd+1)
x2=s'E'
endif
endif
t=t+1
x=x+dx
s=s+dsx
endwhile

x=x0
y=yb
s=LatA
t=1
while (t<nd+2)
if (s<'0')
sa=-s
'draw string 'x' 'y' `1'sa'S'
if (t=1)
y1=sa'S'
endif
if (t=nd+1)
y2=sa'S'
endif
else
if (s!='0')
'draw string 'x' 'y' `1's'N'
else
'draw string 'x' 'y' `1's
endif
if (t=1)
y1=s'N'
endif
if (t=nd+1)
y2=s'N'
endif
endif
t=t+1
x=x+dx
s=s+dsy
endwhile

if (pn=0) ; zca=z0a ; zcb=z0b ; endif
if (pn=1) ; zca=z1a ; zcb=z1b ; endif
if (pn=2) ; zca=z2a ; zcb=z2b ; endif
if (pn=3) ; zca=z3a ; zcb=z3b ; endif
if (pn=4) ; zca=z4a ; zcb=z4b ; endif
if (etc='ett') ; zca=z0a ; zcb=z4b ; endif
ha=gethn(zca)
hb=gethn(zcb)
za=zca-1
zb=zcb-1

if (etc='ett')
  subtit='All Classes : H'za'='ha' to H'zb'='hb
else
  subtit=vn' : Class 'pn' : H'za'='ha' to H'zb'='hb
endif
'draw title 'tit'\`1'subtit

figname='Energy_'dp'_Cross_Sec_'vn'_Class'pn'_'datego'_'pt'_hPa_'caso'_'_Trunc'.png'
say ' '
say ' fig = 'figname
say ' '
fig=_dirfig''figname
'printim 'fig' png white'

return

function gethn(zz)
'set dfile 2'
'set x 1'
'set y 1'
'set z 1'
'display dn.2(x='zz',t=1)'
hh=subwrd(result,4)
if (hh>=1000)
fmt='%6.1f'
endif
if (hh>=100 & hh<1000)
fmt='%6.2f'
endif
if (hh>=10 & hh<100)
fmt='%6.3f'
endif
if (hh>=1 & hh<10)
fmt='%6.4f'
endif
if (hh<1)
fmt='%6.4f'
endif
hn=math_format(fmt,hh)
'set dfile 1'
return hn

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