*  
* Compute Empirical Orthogonal Functions (EOFs)
*
* The current GrADS dimension environment is used for the 
* estimation of the covariance matrix.
* EOFs and PCs are written to separate grads file pairs.
* 
*
* USAGE: [run] eof [options ] <fld>
*
* with 
*    fld:            Input field (GrADS expression).
* Options: 
*   -neof <# eof>  :  Number of EOFs (default: 12) 
*   -npc <# pc>    :  Number of PCs (default: # eof)
*   -tinc <tinc>   :  Time increment (default: 1)
*   -nper <nper>   :  Minimal percent of defined data values required 
*                    for a grid point to be considered in the EOF calculation
*                    (default: 70 %)
*   -o <prefix>    :  Prefix for output files: (default: first 6 characters of <fld>)
*
* Normalization: Variance(PCs)=1
*
* NOTES: 
*        i) The data are not weighted. Use cos(lat*3.141/180) for a simple area weight.
*       ii) The program is not well tested for inputs with data gaps. 
*           Please let me know if you get wrong results with for such input. 
*      iii) Try avoiding GrADS's cyclic continuation for global data. 
*           Otherwise the cyclic points are weighted twice.
*
* (c) Matthias Munnich (munnich@atmos.ucla.edu)
****************************************************
* use MATLAB multivariate EOF by (senya@atmos.umd.edu)
*
* uses WRITECTL GET_UDEF GXSTATE1
* all variables should be regridded to the same grid
* current file grid must correspond to that grid
*****************************************************



function eofnew(args)

'disable fwrite'
*
* get options
*
rc=parsearg(args)
if (rc = 1)
  return
endif

****************
_undef=get_undef()
****************
**********************
printopt()
**********************



*save dimensions *****
dimensions=gxstate1()
**********************

********************************************************
*   WRITE EOF and PC CTLs
********************************************************

**********************
*write EOF ctl
**********************
i=1
str=''
while (i<=_nvars)
str=str' '_fld.i'>EOF'i
i=i+1
endwhile


'set t 1 '_neof
args=str' -o '_outb'_eof 0'
result=writectl(args)
nx=subwrd(result,1)
ny=subwrd(result,2)
nz=subwrd(result,3)
say '*******************************************'

*  restore dimensions ***
gxstate1(dimensions)
*************************

*************************
*write PC ctl
*************************

*say 'tinc='_tinc
'set x 1'
'set y 1'
'set z 1 '_neof
args='lat>PC -o '_outb'_pc 0 '_tinc
result=writectl(args)
nt=subwrd(result,4)
say '*******************************************'

*  restore dimensions ****
gxstate1(dimensions)
**************************

********************************************************
*   ENF OF: WRITE EOF and PC CTLs
********************************************************

********************************************************
* write input data files
********************************************************
'set gxout fwrite'
'set fwrite eof_in.gad'

********************************************************
*   set undef: default is -9.99e8 in v2
********************************************************
'q config'
res=sublin(result,1)
res1=subwrd(res,2)
ver=substr(res1,2,1)
if ver=2
'set undef '_undef
endif

******* adjust to integer grid points ******************
xl=subwrd(dimensions,1)
xr=subwrd(dimensions,2)
yb=subwrd(dimensions,5)
yt=subwrd(dimensions,6)
zb=subwrd(dimensions,9)
zt=subwrd(dimensions,10)
iz1=subwrd(dimensions,11)
iz2=subwrd(dimensions,12)
t1=subwrd(dimensions,15)
t2=subwrd(dimensions,16)

*adjust range to grid 
'set lon 'xl
xl1=subwrd(result,4)
'set lon 'xr
xr1=subwrd(result,4)

'set lat 'yb
yb1=subwrd(result,4)
'set lat 'yt
yt1=subwrd(result,4)

'set lev 'zb
zb1=subwrd(result,4)
'set lev 'zt
zt1=subwrd(result,4)

'set lon 'xl1' 'xr1
'set lat 'yb1' 'yt1
'set lev 'zb1' 'zt1
'set t 't1' 't2
**** end adjust to grid ****************

say '      Writing data to transfer file...'

j=1
while (j<=_nvars)
tt=t1
tdim=0
while (tt<=t2)
   'set t 'tt
   'd '_fld.j
   tt=tt+_tinc
   tdim=tdim+1
endwhile
j=j+1
endwhile
'disable fwrite'
say '      'tdim ' time steps written'

'set gxout contour'
say '*******************************************'

********************************************************
* END OF: write input data files
********************************************************

*  restore dimensions ****
gxstate1(dimensions)
**************************
* write info file
dum=write('eof_in.inf','% nx ny nz nt nvars neof nper undef file_name')
dum=write('eof_in.inf',nx' 'ny' 'nz' 'nt' '_nvars' '_neof' '_nper' '_undef' '_outb, append)
dum=close('eof_in.inf')
**************************

************************************************************
say '        *** computing EOFs ***'
************************************************************
'!matlab -nojvm -r "eofgrads; quit"'
*'run_matlab_eof'
************************************************************

****** cleaning files **************************************
'!rm eof_in.gad eof_in.inf'
************************************************************
return 
****************************************************

***************************
*****************************************************
function parsearg(arg)
* Defaults
_neof=12
_nper=70
_npc=-1
_verb=0
_tinc=1
_outb=''

i=1
word=subwrd(arg,i)
if(word = '')
  usage()
  return 1
endif

while (substr(word,1,1) = '-') 
  if (word = '-neof')
     i=i+1
     _neof=subwrd(arg,i)
  else
    if (word = '-npc') 
       i=i+1
       _npc=subwrd(arg,i)
    else
      if (word = '-nper')
        i=i+1
        _nper=subwrd(arg,i)
      else
        if (word = '-o' | word = '-outname')
          i=i+1
          _outb=subwrd(arg,i)
        else
          if(word = '-v' | word = '-verbos')
            _verb = 1
          else
            if(word = '-tinc')
              i=i+1
              _tinc=subwrd(arg,i)
            else
              say ' Unknown option: "'word'"'
              usage()
              return 1
            endif
          endif              
        endif
      endif
    endif
  endif
  i=i+1
  word=subwrd(arg,i)
endwhile

j=1
fld=subwrd(arg,i)
if(_outb= '') 
   _outb=substr(fld,1,6)
endif
if(fld = '')
  say ' Fatal: Missing field expression'
  usage()
  return 1
endif

while (fld!='')
_fld.j=fld
i=i+1
j=j+1
fld=subwrd(arg,i)
endwhile
_nvars=j-1


if(_verb>0)
  printopt()
endif
return 0
*******************************************************
*****************************************************
function printopt()
say '***************************************'
say '      Number of EOFs:         ' _neof
*say '      Number of PCs:          ' _npc
say '      Time increment:         ' _tinc
say '      Required % of good data:' _nper
say '      Output files basename:  ' _outb
say '      Field undef:            ' _undef
*say '      Field file is XDF:      ' _xdf
say '      Number of fields:       ' _nvars
ivar=1
while (ivar<=_nvars) 
say '      EOF'ivar' of :               ' _fld.ivar
ivar=ivar+1
endwhile
* say '      Verbose:' _verb
say '***************************************'
return
*****************************************************

function usage()
say ' EOF: Compute Empirical Orthogonal Functions (EOFs)'
say '                eof.gs Version 0.14'
say ' '
say ' USAGE: eof [options] fld'
say ' '
say ' with '
say '    fld:            Input field (GrADS expression).'
say ' '
say ' OPTIONS: '
say '   -neof <# eof> :  Number of EOFs (default: 12)'
say '   -npc <# pc>   :  Number of PCs (default: # eof)'
say '   -tinc <tinc>  :  Time increment (default: 1)'
say '   -nper <nper>  :  Minimal percent of defined data values required '
say '                    for a grid point to be considered in the EOF calculation'
say '                    (default: 70 %)'
say '   -o <prefix>   :  Prefix for files: (default: first 6 characters of <fld>)'
say ' NOTES: '
say '        i) The data are not weighted. Use cos(lat*3.141/180) for a simple area weight.'
say '       ii) The program is not well tested for inputs with data gaps.'
say '           Please let me know if you get wrong results with for such input.'
say ' (c) Matthias Munnich (munnich@atmos.ucla.edu)'
return
*****************************************************

function get_undef()
'q ctlinfo'
i=1
while i<=100000
res=sublin(result,i)
qq=subwrd(res,1)
if (qq='undef' | qq='UNDEF' |qq='Undef')
undef=subwrd(res,2)
return undef
endif
i=i+1
endwhile
return
******************************************************

******************************************************
function gxstate1(args)
*
*  Get or restore dimensions 
* get dimension USAGE: gxstate1
* restore dimension USAGE: gxstate1 dimensions
* dimeansions is the output of get dimension by gxstate1()
******************
if args='args' | args='' 

'query dim'
dinf = result
lx = sublin(dinf,2)
ly = sublin(dinf,3)
lz = sublin(dinf,4)
lt = sublin(dinf,5)
if ( subwrd(lx,7) = 'to')
  lons = subwrd(lx,6)
  lone = subwrd(lx,8)
  xs = subwrd(lx,11)
  xe = subwrd(lx,13)
else
  lons = subwrd(lx,6)
  lone = subwrd(lx,6)
  xs = subwrd(lx,9)
  xe = subwrd(lx,9)
endif
if ( subwrd(ly,7) = 'to')
  lats = subwrd(ly,6)
  late = subwrd(ly,8)
  ys = subwrd(ly,11)
  ye = subwrd(ly,13)
else
  lats = subwrd(ly,6)
  late = subwrd(ly,6)
  ys = subwrd(ly,9)
  ye = subwrd(ly,9)
endif
if ( subwrd(lz,7) = 'to')
  levs = subwrd(lz,6)
  leve = subwrd(lz,8)
  zs = subwrd(lz,11)
  ze = subwrd(lz,13)
else
  levs = subwrd(lz,6)
  leve = subwrd(lz,6)
  zs = subwrd(lz,9)
  ze = subwrd(lz,9)
endif
if ( subwrd(lt,7) = 'to')
  tims = subwrd(lt,6)
  time = subwrd(lt,8)
  ts = subwrd(lt,11)
  te = subwrd(lt,13)
else
  tims = subwrd(lt,6)
  time = subwrd(lt,6)
  ts = subwrd(lt,9)
  te = subwrd(lt,9)
endif

return(lons' 'lone' 'xs' 'xe' 'lats' 'late' 'ys' 'ye' 'levs' 'leve' 'zs' 'ze' 'tims' 'time' 'ts' 'te)

********end get dimension***********
else
'set lon 'subwrd(args,1)' 'subwrd(args,2)
'set lat 'subwrd(args,5)' 'subwrd(args,6)
'set lev 'subwrd(args,9)' 'subwrd(args,10)
*'set time 'subwrd(args,13)' 'subwrd(args,14)
'set t 'subwrd(args,15)' 'subwrd(args,16)
endif
*****************************************************

function writectl(args)

if args='args' | args=''
say 'set area and time'
say 'USAGE: writevar var1[>newname1] var2[>newname2] ...'  
say '...-o fname(no extension) [idat] [tinc]'
say '[...] optional, default is the same name, idat=1, tinc=1'
return
endif

*read list variables
ans='qq'
i=1
while (ans!='' & ans!='-o')
var2.i=subwrd(args,i)
i=i+1
ans=subwrd(args,i)
endwhile
nvars=i-1

*read output file name or set default
if (ans='')
fn='qq#'
else
fn=subwrd(args,i+1)
endif 

*idat=0 only CTL, idat=1 CTL and DATA
idat=subwrd(args,i+2)
if idat=''
idat=1
endif

*tinc
tinc=subwrd(args,i+3)
if tinc=''
tinc=1
endif

say 'OUT_FILE: 'fn

*check for variable rename**********
i=1
while(i<=nvars)
vl=strlen(var2.i)
j=1
i1=vl

while j<=vl
if substr(var2.i,j,1)='>'
i1=j
endif

j=j+1
endwhile

if i1=vl
var.i=var2.i
varn.i=var2.i
else
var.i=substr(var2.i,1,i1-1)
varn.i=substr(var2.i,i1+1,vl-i1)
endif
say var.i '-->'varn.i
i=i+1
endwhile



*save settings
dimensions=gxstate1()

xl=subwrd(dimensions,1)
xr=subwrd(dimensions,2)
yb=subwrd(dimensions,5)
yt=subwrd(dimensions,6)
zb=subwrd(dimensions,9)
zt=subwrd(dimensions,10)
iz1=subwrd(dimensions,11)
iz2=subwrd(dimensions,12)
t1=subwrd(dimensions,15)
t2=subwrd(dimensions,16)

*adjust range to grid 
'set lon 'xl
xl1=subwrd(result,4)
'set lon 'xr
xr1=subwrd(result,4)

'set lat 'yb
yb1=subwrd(result,4)
'set lat 'yt
yt1=subwrd(result,4)

'set lev 'zb
zb1=subwrd(result,4)
'set lev 'zt
zt1=subwrd(result,4)

'set lon 'xl1' 'xr1
'set lat 'yb1' 'yt1
'set lev 'zb1' 'zt1
'set t 't1' 't2

*extract CTL file information
'q dims'
xline=sublin(result,2)
yline=sublin(result,3)
zline=sublin(result,4)
tline=sublin(result,5)
xstatus=subwrd(xline,3)
ystatus=subwrd(yline,3)
zstatus=subwrd(zline,3)
tstatus=subwrd(tline,3)

if xstatus='fixed'
xset='XDEF 1 LEVELS '%subwrd(xline,6)
endif
if ystatus='fixed'
yset='YDEF 1 LEVELS '%subwrd(yline,6)
endif
if zstatus='fixed'
zset='ZDEF 1 LEVELS '%subwrd(zline,6)
levsize=1
endif
if tstatus='fixed'
tset='TDEF 1 LINEAR '%subwrd(tline,6)
endif

 

'set gxout stat'

if xstatus!='fixed'
'set lon 'xl1' 'xr1
'set y 1'
'set z 1'
'set t 1'


'd ' var.1
sline=sublin(result,5)



xline=sublin(result,3)
xsize=subwrd(sline,3)
xscale=subwrd(xline,7)


lonsize=xsize
if(xscale = 'Linear'); xsize=2; endif;
i=8
while(i<xsize+9)
  j=i-7
  xg.j=subwrd(xline,i)
  i=i+1
endwhile
endif


if ystatus!='fixed'
'set lat 'yb1' 'yt1
'set x 1'
'set z 1'
'set t 1'

'd ' var.1
sline=sublin(result,5)
yline=sublin(result,3)
ysize=subwrd(sline,3)
yscale=subwrd(yline,7)


latsize=ysize
if(yscale = 'Linear'); ysize=2; endif;
i=8
while(i<ysize+9)
  j=i-7
  yg.j=subwrd(yline,i)
  i=i+1
endwhile
endif


if zstatus!='fixed'
'set lev 'zb1' 'zt1
'set y 1'
'set x 1'
'set t 1'

'd ' var.1
sline=sublin(result,5)
zline=sublin(result,3)
zsize=subwrd(sline,3)
zscale=subwrd(zline,7)


levsize=zsize
if(zscale = 'Linear'); zsize=2; endif;
i=8
while(i<zsize+9)
  j=i-7
  zg.j=subwrd(zline,i)
  i=i+1
endwhile
endif



'set x 1'
'set z 1'
'set y 1'
if tstatus='fixed'
t22=t2+1
else
t22=t2
endif
'set t 't1' 't22

*say var.1' 't1' 't22



*'q dims'
*say result
'buf='var.1
'd  buf'
***wrong if variable is an expression of variables with different undefs****
*uline=sublin(result,6)
*undef=subwrd(uline,4)

*say result


tline=sublin(result,3)
sline=sublin(result,5)
tsize=subwrd(sline,3)
tscale=subwrd(tline,7)


timesize=tsize

if tstatus='fixed'
timesize=1
endif


if(tscale = 'Linear'); tsize=2; endif;
*say tline
i=8
while(i<tsize+9)
  j=i-7
  tg.j=subwrd(tline,i)
*  say 'tg.'j' 'tg.j
  i=i+1
endwhile

*write CTL
*********************
fn1=get_filename(fn)
*********************
dum=write(fn'.ctl','DSET ^'fn1'.dat')
dum=write(fn'.ctl','UNDEF '_undef, append)
dum=write(fn'.ctl','TITLE 'fn, append)
*X-LINE
if (xstatus='fixed')
dum=write(fn'.ctl',xset,append)
nx=1
else
if (xscale='Linear')
dum=write(fn'.ctl','XDEF      'lonsize' 'xscale' 'xg.1' 'xg.2,append)
nx=lonsize
else
dum=write(fn'.ctl','XDEF      'lonsize' 'xscale,append)
nx=lonsize
  n5=xsize/5
  i=1
  while(i<xsize+1)
    a1=xg.i
    j=i+1; a2=xg.j
    j=i+2; a3=xg.j
    j=i+3; a4=xg.j
    j=i+4; a5=xg.j
    dum=write(fn'.ctl',a1' 'a2' 'a3' 'a4' 'a5,append)
    i=i+5
  endwhile
  a1=''
  a2=''
  a3=''
  a4=''
  a5=''
  endif
endif  
  
*Y-LINE
if ystatus='fixed'
dum=write(fn'.ctl',yset,append)
ny=1
else
if (yscale='Linear')
dum=write(fn'.ctl','YDEF      'latsize' 'yscale' 'yg.1' 'yg.2,append)
ny=latsize
else
dum=write(fn'.ctl','YDEF      'latsize ' 'yscale,append)
ny=latsize
  n5=ysize/5
  i=1
  while(i<ysize+1)
    a1=yg.i
    j=i+1; a2=yg.j
    j=i+2; a3=yg.j
    j=i+3; a4=yg.j
    j=i+4; a5=yg.j
    dum=write(fn'.ctl',a1' 'a2' 'a3' 'a4' 'a5,append)
    i=i+5
  endwhile
  a1=''
  a2=''
  a3=''
  a4=''
  a5=''
endif
endif
  
*Z-LINE
if zstatus='fixed'
dum=write(fn'.ctl',zset,append)
nz=1
else
if (zscale='Linear')
dum=write(fn'.ctl','ZDEF      'levsize' 'zscale' 'zg.1' 'zg.2,append)
nz=levsize
else
dum=write(fn'.ctl','ZDEF      'levsize ' 'zscale,append)
nz=levsize
  n5=zsize/5
  i=1
  while(i<zsize+1)
    a1=zg.i
    j=i+1; a2=zg.j
    j=i+2; a3=zg.j
    j=i+3; a4=zg.j
    j=i+4; a5=zg.j
    dum=write(fn'.ctl',a1' 'a2' 'a3' 'a4' 'a5,append)
    i=i+5
  endwhile
  a1=''
  a2=''
  a3=''
  a4=''
  a5=''  
endif
endif

*T-LINE
if tstatus='fixed'
dum=write(fn'.ctl',tset' 'tg.2,append)
nt=1
else
if tinc=1
dum=write(fn'.ctl','TDEF      'timesize' 'tscale' 'tg.1' 'tg.2,append)
nt=timesize
else
ntime=math_int((timesize-1)/tinc)+1
tstrl=strlen(tg.2)
tint=tinc%substr(tg.2,tstrl-1,2)
dum=write(fn'.ctl','TDEF      'ntime' 'tscale' 'tg.1' 'tint,append)
nt=ntime
endif
endif


dum=write(fn'.ctl','VARS 'nvars,append)
i=1
while(i<nvars+1)
dum=write(fn'.ctl',varn.i' 'levsize' 99 **'var.i,append)
i=i+1
endwhile

dum=write(fn'.ctl','ENDVARS ',append)

if idat=1
*write data
'set lon 'xl1' 'xr1
'set lat 'yb1' 'yt1
*'set lev 'zb1' 'zt1
'set z 'iz1

'set gxout fwrite'
'set fwrite 'fn'.dat'
t=t1
while(t<=t2)
'set t 't
*say t
i=1
while (i<=nvars)
iz=iz1
while(iz<=iz2)
'set z 'iz
'd 'var.i
iz=iz+1
endwhile
i=i+1
endwhile
t=t+tinc
endwhile
'disable fwrite'
'set gxout contour'

endif
gxstate1(dimensions)
return (nx' 'ny' 'nz' 'nt' ')
**** end write var *****************************************
function get_filename (args)
if args=''
say 'USAGE: get_filename path/filename '
return
endif

fn=subwrd(args,1)
n=strlen(fn)


path=get_path(fn)
if path=''
n1=0
else
n1=strlen(path)
endif



filename=substr(fn,n1+1,n-n1)
return (filename)
*************************************
function get_path (args)
if args=''
say 'USAGE: get_path filename '
return
endif

fn=subwrd(args,1)

qq='a'
i2=0
i=1
while (qq!='')
qq=substr(fn,i,1)
if(qq='/' | qq='\')
i2=i
endif
i=i+1
endwhile

if i2!=0
pathname=substr(fn,1,i2)
else
pathname=''
endif
return(pathname)
**************************************