function ref (args)

Mmax=subwrd(args,1)
if (Mmax='-h')
  say ' '
  say ' Plots the Reinforcement Points of Energy for the Case of Barotropic'
  say ' Mixed Waves 3 and 1 which gives the Best and Suitable Graphic'
  say ' '
  say ' Usage:'
  say ' Reinforcement [Mmax]'
  say ' '
  say ' Default: Mmax = 126'
  say ' '
  say ' Must be runned under grads landscape'
  say ' '
  say ' arg meaning:'
  say ' '
  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 ' '
  say '  OBS:'
  say '      Mixed Rossby-Gravity: mx : l=Mmax+17, isi=2, ln=0'
  say '      Wave Numbers for Reinforcement: s1=3 and s2=1'
  say ' '
  return
endif
if (Mmax='') ; Mmax=126 ; endif

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

* Getting Kmax (number of pressure levels) and referece case csvm
'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)
'!rm -f 'filecsdt

* 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)
if (Kmax < 100)
  Trunc=Rhomb'L0'Kmax
else
  Trunc=Rhomb'L'Kmax
endif

say ' '
say ' Mmax = 'Mmax'   Rhomb = 'Rhomb'   Trunc = 'Trunc
say ' '

cs='mx'
sn1='s03'
sn2='s01'
ln='l00'
hex='h001'

'reinit'
'set display grey white'
'clear'

refname='Reinforcement_'cs'_'sn1'_'sn2'_'ln'_'Trunc'_'csvm'_'hex
dgs='./dgs/'refname'.dgs'
say ' '
say ' 'refname
say ' '

rec=read(dgs)
lrec=sublin(rec,2)
ctit=lrec

rec=read(dgs)
lrec=sublin(rec,2)
eps=subwrd(lrec,1)
hn=subwrd(lrec,2)

rec=read(dgs)
lrec=sublin(rec,2)
s1=subwrd(lrec,1)
n1=subwrd(lrec,2)

rec=read(dgs)
lrec=sublin(rec,2)
s2=subwrd(lrec,1)
n2=subwrd(lrec,2)

rec=read(dgs)
lrec=sublin(rec,2)
to=subwrd(lrec,1)
tm=subwrd(lrec,2)
lo=subwrd(lrec,3)
lm=subwrd(lrec,4)

rec=read(dgs)
lrec=sublin(rec,2)
dt=subwrd(lrec,1)
dl=subwrd(lrec,2)
cg=subwrd(lrec,3)

xo=2.0
yo=1.0-0.4
xm=10.0
ym=6.0-0.4
im=8
dx=(xm-xo)/im
jm=8
dy=(ym-yo)/jm
tk=0.2
ss=0.2
st=0.85*ss
yt=8.3
sb=0.75*ss
sx=2.25*st
sy=sb
mrksiz=0.15
ax=(xm-xo)/(lm-lo)
bx=xo-ax*lo
ay=(ym-yo)/(tm-to)
by=yo-ay*to

'set line 1 1 6'
'draw line 'xo' 'yo' 'xm' 'yo
'draw line 'xo' 'ym' 'xm' 'ym
'draw line 'xo' 'yo' 'xo' 'ym
'draw line 'xm' 'yo' 'xm' 'ym

'set strsiz 'ss
'set string 1 c 6 0'
'draw string 'xo' 'yo-ss' -`3p'
'draw string 'xo+2*dx' 'yo-ss' -`3p`0/2'
'draw string 'xo+4*dx' 'yo-ss' 0'
'draw string 'xo+6*dx' 'yo-ss' `3p`0/2'
'draw string 'xm' 'yo-ss' `3p'
'draw string 'xo+4*dx' 'yo-2.25*ss' `3l'
i=1
while (i <= im-1)
'draw line 'xo+i*dx' 'yo' 'xo+i*dx' 'yo+tk
'draw line 'xo+i*dx' 'ym' 'xo+i*dx' 'ym-tk
i=i+1
endwhile

'set string 1 r 6 0'
'draw string 'xo-ss' 'yo' 0'
dj=tm/jm
j=1
while (j <= jm-1)
'draw line 'xo' 'yo+j*dy' 'xo+tk' 'yo+j*dy
'draw line 'xm' 'yo+j*dy' 'xm-tk' 'yo+j*dy
'draw string 'xo-ss' 'yo+j*dy' 'dj*j
j=j+1
endwhile
'draw string 'xo-ss' 'ym' 'tm
'set string 1 c 6 90'
'draw string 'xo-5*ss' 'yo+4*dy' 2 `3Z`0 t'

mm=100
while (mm!=-9999)
'set line 1 2 6'
rec=read(dgs)
lrec=sublin(rec,2)
mm=subwrd(lrec,1)
if (mm!=-9999)
l1=subwrd(lrec,2)
x1=ax*l1+bx
t1=subwrd(lrec,3)
y1=ay*t1+by
rec=read(dgs)
lrec=sublin(rec,2)
mm=subwrd(lrec,1)
l2=subwrd(lrec,2)
x2=ax*l2+bx
t2=subwrd(lrec,3)
y2=ay*t2+by
'draw line 'x1' 'y1' 'x2' 'y2
rec=alfa(x1,y1,x2,y2,'m',mm,sx,sy,sb)
endif
endwhile
rec=read(dgs)
lrec=sublin(rec,2)
m1tit=subwrd(lrec,1)

mm=100
while (mm!=-9999)
'set line 1 6 6'
rec=read(dgs)
lrec=sublin(rec,2)
mm=subwrd(lrec,1)
if (mm!=-9999)
l1=subwrd(lrec,2)
x1=ax*l1+bx
t1=subwrd(lrec,3)
y1=ay*t1+by
rec=read(dgs)
lrec=sublin(rec,2)
mm=subwrd(lrec,1)
l2=subwrd(lrec,2)
x2=ax*l2+bx
t2=subwrd(lrec,3)
y2=ay*t2+by
'draw line 'x1' 'y1' 'x2' 'y2
rec=alfa(x1,y1,x2,y2,'m',mm,sx,sy,sb)
endif
endwhile
rec=read(dgs)
lrec=sublin(rec,2)
m2tit=subwrd(lrec,1)

mm=100
while (mm!=-9999)
'set line 1 1 6'
rec=read(dgs)
lrec=sublin(rec,2)
mm=subwrd(lrec,1)
if (mm!=-9999)
l1=subwrd(lrec,2)
x1=ax*l1+bx
t1=subwrd(lrec,3)
y1=ay*t1+by
rec=read(dgs)
lrec=sublin(rec,2)
mm=subwrd(lrec,1)
l2=subwrd(lrec,2)
x2=ax*l2+bx
t2=subwrd(lrec,3)
y2=ay*t2+by
'draw line 'x1' 'y1' 'x2' 'y2
rec=alfa(x1,y1,x2,y2,'j',mm,sx,sy,sb)
endif
endwhile
rec=read(dgs)
lrec=sublin(rec,2)
jtit=subwrd(lrec,1)

'set strsiz 'ss
'set line 1 5 6'
l1=0
x1=ax*l1+bx
t1=0
y1=ay*t1+by
l2=0
x2=ax*l2+bx
t2=dt
y2=ay*t2+by
'draw line 'x1' 'y1' 'x2' 'y2
xs=(x1+x2)/2+2*ss
ys=(y1+y2)/2+ss
'set string 1 c 6 0'
'draw string 'xs' 'ys' `3D`0t`b1,2'
l1=0
x1=ax*l1+bx
t1=dt
y1=ay*t1+by
l2=dl
x2=ax*l2+bx
t2=dt
y2=ay*t2+by
'draw line 'x1' 'y1' 'x2' 'y2
xs=(x1+x2)/2
ys=(y1+y2)/2-ss
'set string 1 c 6 0'
'draw string 'xs' 'ys' `3Dl`0`b1,2'

jj=100
while (jj!=-9999)
rec=read(dgs)
lrec=sublin(rec,2)
jj=subwrd(lrec,1)
mm=subwrd(lrec,2)
if (jj!=-9999)
ll=subwrd(lrec,3)
xx=ax*ll+bx
tt=subwrd(lrec,4)
yy=ay*tt+by
'draw mark 3 'xx' 'yy' 'mrksiz
endif
endwhile

'set strsiz 'st
'set string 1 l 6 0'
xs=xo/2+st
ys=yt
'draw string 'xs' 'ys' Reinforcement Points (  ) between the Waves 'ctit' - 'csvm':'
xx=xs+19.5*st
yy=ys
'draw mark 3 'xx' 'yy' 'mrksiz
'set string 1 c 6 0'
xs=xo+(xm-xo)/2-xo/4
ys=ys-2*st
'draw string 'xs' 'ys' s`b1`n='s1', `3w`0`b1`n='n1' Wave I, s`b2`n='s2', `3w`0`b2`n='n2' Wave II, `3e`0`bn`n='eps' H`bn`n='hn
ys=ys-2*st
'draw string 'xs' 'ys' Cg = (`3w`0`b2`n-`3w`0`b1`n)/(s`b2`n-s`b1`n) = `3Dl`0`b1,2`n/`3D`0t`b1,2`n = 'cg
ys=ys-2*st
'draw string 'xs' 'ys' `3Dl`0`b1,2`n = 2`3p`0(`3w`0`b2`n-`3w`0`b1`n)/(s`b1`n`3w`0`b2`n-s`b2`n`3w`0`b1`n) = 'dl
ys=ys-2*st
'draw string 'xs' 'ys' `3D`0t`b1,2`n = 2`3p`0(s`b2`n-s`b1`n)/(s`b1`n`3w`0`b2`n-s`b2`n`3w`0`b1`n) = 'dt
'set string 1 c 6 0'
ys=ys-2*st
'draw string 'xs' 'ys' Continuous Lines: Contant Phase Diferences, Waves I and II, 'jtit
ys=ys-2*st
'draw string 'xs' 'ys' Dashed Lines: Wave I Crest Positions, 'm1tit
ys=ys-2*st
'draw string 'xs' 'ys' Dashed-Doted Lines: Wave II Crest Positions, 'm2tit

'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-0.4
ys=yb-1.3

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='../png_ERA_5/'refname'.png'
say ' '
say ' 'fig
say ' '
'printim 'fig' png white'

return

function alfa(x1,y1,x2,y2,h,m,sx,sy,sb)

if (x1=x2) ; return ; endif 
a=(y2-y1)/(x2-x1)
b=y1-a*x1
xm=x1+3.25*(x2-x1)/5
ym=a*xm+b
xa=xm-sx
xb=xm+sx
ya=ym-sy
yb=ym+sy
'set line 0 1 6'
'draw recf 'xa' 'ya' 'xb' 'yb
'set line 1 1 6'
'draw rec 'xa' 'ya' 'xb' 'yb
'set strsiz 'sb
'set string 1 c 6 0'
'draw string 'xm' 'ym' 'h'='m

return

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