#!/bin/sh
#example calculation of S/N EOFs
set -x

ingrid<<EOF
\beginingrid

% definition to change  the names of grids
/renamegrids {time /T renameGRID lat /Y renameGRID lon /X renameGRID} def

% mask out land
/ccm_mask {
SOURCES .IRI .MP .RESEARCH .ATMOSPHERE .GLOBAL .CCM .mem01 .ORO
  0 maskgt
  [time]average 1 add 
} def

% compute ensemble mean

SOURCES .IRI .MP .RESEARCH .ATMOSPHERE .GLOBAL .CCM
     a: .mem01 .TAUX
    :a: .mem02 .TAUX add
    :a: .mem03 .TAUX add
    :a: .mem04 .TAUX add
    :a
4 div
ccm_mask mul
/name /TAUX def
renamegrids
(ensemble_mean.cdf) writeCDF	% write full field

SOURCES .IRI .MP .RESEARCH .ATMOSPHERE .GLOBAL .CCM
     a: .mem01 .TAUX
    :a: .mem02 .TAUX add
    :a: .mem03 .TAUX add
    :a: .mem04 .TAUX add
    :a
4 div
ccm_mask mul
/name /TAUX def
renamegrids
yearly-anomalies
(ensemble_meana.cdf) writeCDF	% write anomalies

EOF


# compute EOFs of ensemble mean anamolies--not needed for S/N calculation.
# nice to have
ingrid<<EOF &
\beginingrid
(ensemble_meana.cdf) readfile
[X Y] [T] svd
(mean_eofs.cdf) writeCDF
EOF


# compute deviations from the ensemble mean
# ugly!
ingrid<<EOF
\beginingrid
/renamegrids {time /T renameGRID lat /Y renameGRID lon /X renameGRID} def
/ensmean { (ensemble_mean.cdf) readfile } def
concatstreams: 
  SOURCES .IRI .MP .RESEARCH .ATMOSPHERE .GLOBAL .CCM 
   a: .mem01 .TAUX renamegrids ensmean .TAUX sub
  nrdim SetStreamIndex* X Y /M unitless 0 1 NewGRID T *
  :a: .mem02 .TAUX renamegrids ensmean .TAUX sub
  nrdim SetStreamIndex* X Y /M unitless 0 2 NewGRID T *
  :a: .mem03 .TAUX renamegrids ensmean .TAUX sub
  nrdim SetStreamIndex* X Y /M unitless 0 3 NewGRID T *
  :a: .mem04 .TAUX renamegrids ensmean .TAUX sub
  nrdim SetStreamIndex* X Y /M unitless 0 4 NewGRID T *
  :a
:concatstreams
(deviations.cdf) writeCDF

EOF

#compute eofs of deviations (noise eofs)
ingrid<<EOF 
\beginingrid
(deviations.cdf) readfile
[X Y] [T M] svd
(deviation_eofs.cdf) writeCDF
EOF

# compute "filter"
ingrid<<EOF
\beginingrid
% factor to normalize eofs to unity
/fac 64. 128. mul sqrt def
(deviation_eofs.cdf) readfile 
  a: .Ss  fac div :a: .sv fac mul :a
   div 
(ensemble_meana.cdf) readfile .TAUX
T /time renameGRID
  mul   [X Y] sum
/name /filter def
ev /J renameGRID
time /T renameGRID
(filter.cdf) writeCDF
EOF

# compute S/N eofs
# Jmax = # of noise eofs to use
Jmax=100
ingrid<<EOF &
\beginingrid
ds:
(filter.cdf) readfile .filter
J 1 $Jmax RANGE
[J] [T] svd 
a: .sv   $Jmax sqrt  mul
% factor to fix sv normalization
  /comment_sn_eof ($Jmax noise EOFs used.) def
:a: .Ts :a 
/name /Ts def
/comment_sn_eof ($Jmax noise EOFs used.) def
$Jmax sqrt mul
dup
CopyStream
dup dup mul [T] sum div
(ensemble_meana.cdf) readfile .TAUX
mul [T] sum
/name /Ss def
/long_name /structures def
/comment_sn_eof ($Jmax noise EOFs used.) def
:ds
/comment_sn_eof ($Jmax noise EOFs used.) def
(s2n_eofs.cdf) writeCDF
EOF

#compute filtered signal
# Smax = # of S/N eofs to use to reconstruct the signal

Smax=10
ingrid<<EOF 
\beginingrid
(s2n_eofs.cdf) readfile
  ev 1 $Smax RANGE
  a: .Ss     :a: .Ts :a
  mul
  [ev]sum
2 RECHUNK 
/name /TAUX def
/long_name /TAUX def
/comment_signal ($Smax S/N EOFs used.) def
(signal.cdf) writeCDF
EOF