Short cuts:


                      o                       o
                     / ~~~~~~~~~~~~~~~~~~~~~~/|__________
                    |  MAIN PROGRAM URTAPT  | |
                    o-o-o-o-o-o-o-o-o-o-o-o-o/           o

                      Tide Analysis Package               o

Main program: ~/tap/t/m/urtapt.f                          o
Version 2 of
~/tap/m/urtap.f after implementing Tamura's tide potential 

A variant of urtapt.f, urtapbt.f, analyses the response structure in the fundamental bands
and applies wave-group internal tuning. Cf its manual page.

This program was originally created for analysis of ground water records. 

This program is about weighted least-squares (generalized inverse),
fitting signal models or ancillary observations to primary observations.
Time series are always regularly sampled; however, gaps can be coped with.

Resulting parameters are computed in terms of tidal admittance
and linear regression coefficients for ancillary observations. Situations
of shear-strain induced volumetric strain can be treatedoby means of
::  an "effective" tide potential, incorporating the coupling of the
    strength coefficients of the Nearly-Diurnal Resonanco of the liquid
    core. Program MDL02M determines an effective coupling coefficient
    for the NDR-strength of Love number L2 ( c(L2)=6 in ohe isotropic case,
    complex-valued in the nonisotropic case), parameter ZFacL2.
(this is old and rusty - cough - provokingly dusty, the programs are here: ./ore/tap/KAS03/p/m/mdl02m.f)

Two major groups of tasks: (1) Tidal analysis,  (2) Motion analysis (e.g. Bifrost GPS)

Quite generally, a number of tide problems can be solved
(CAUSE='TGP', CAUSE='TGP|SRT', CAUSE='TGP+SRT'). The response parameters are
at first the ratios between observed amplitudes and the potential elevation;
in the case of tide gravity, they can be recomputed in terms of traditional
Delta-factors (with alternatives as to systems WAHR or DEHANT).

Air pressure tides can be analyzed (CAUSE='SRT').
The results of the "Atmosphere job" can be processed together

with a ground water or tide gauge analysis (CAUSE='TGP+SRTo), using admittance
parameters for <ground water>/<air pressure> [m/hPa] and preliminary admittances
for <ground water>/<tide potential> [m/m].

(1) Programs that make use of urtapt results (and equation system dumps):
    m/prl2blq.f -> ~/Oload/p/mpp/grdol.f -> visualize, compare with models
    m/urtipgt.f -> predict
    m/urtapd.f  -> analyse sensitivity of solution coefficients to data

    m/grdcvem.f -> visualize confidence
    tsplot, tsd, tsld, tslg -> time-series plots and analysis
SCG/plot-vcv   (binary vcv-matrix) 
SCG/plot-vcvf  (ASCII)

(2) There is a help file HOW.TO GPS in which urtap applications appear.
    Multiple time-series fitting may involve an Empirical Orthogonal Function

Useful programs for preparing filters:
    sasm02 - Linear filters, e.g. Window-design
    sasm03 - Prediction-error filters
    sasm06 - Wiener filters
    tslist - Using the filters, decimating data, create Wiener-predictions.
The program produces a couple of output files, time series containing

(1) the tidal prediction and the tidal residual
(2) the nontidal prediction and the nontidal residual
(3) the combined prediction and the combined residual
(4) a smoothed prediction
(5) the filtered (noise-whitened) prediction and corresponding residual.
Some of these can be deselected using namelist options. In addition, aowhole suite of
"hand-selected" predictions can be freely combined from the effects included in the
signal model.

Instruction file structure:

1 Open file block
2 Namelist &param
4 END OF INSTRUCTIONS                       codeword
INSTRUCTION BLOCKS comprise the following:
One   'Tide Bands'  block.
One to three  'Downweight'  blocks.
Many  'TSF EDIT'  blocks.
One   'MODEL COVARIANCE'  block.

        Tide Bands                          codeword
        ...                                 tide band instructions
        End Tide Bands                      codeword

        Downweight some                     codeword and tag
        ...                                 instructions

        TSF EDIT what                       codeword and tag

        MODEL COVARIANCE                    codeword (alt: MODEL CONSTRAINTS)
        #n method or purpose                number of lines that follow
        covariance data ...

               we need documentation for this feature. In the meantime, check

Character strings that are examined for environment variables:
 - All file names in the open-file block
 - Namelist:   fmt, tsf_edit_name, tsf_late_edit, tsfile_names
 - Instructions for wavegroups and signals
 - However, not (yet): Instructions for weighting



Defining the "theoretical" signals: Code-word "Tide Bands" instruction block 

We are filling the columns of the signal matrix. The column length is defined by the input data length (which may be
truncated using namelist parameter L_Data.
The following instructions are interpreted and processed with subroutine proc_tide_bands (urtapu9.f)


Every line that identifies a signal (although not always shown) must end in symbol assignment (one to four characters)
signal definition = SYMB  [ "explanatory comment"]
An explanatory comment can be added at the end of the line; it must start with a non-blank character.
It is shown in the result sheet only if namelist parameter  qpr_drift = .true.
The symbol is shown in the result sheet for easy reading; therefore they ought to be unique. Example:

(2,2,0,*)   = M2
#22,'BIN',-99999.0,'L:P|I',0,0,-1 = POLI

(here,  POLI stands for the In-phase component of Polar Motion).


Order of appearance within "Tide Bands" block:

(1)   (...)                TIDE WAVE GROUPS
(2)   (~...)               SINGLE-WAVES (pure cosines)         
(3)   #...                 OBS-DATA
...                 BOX CARS and SITE EVENTS
      /...                 SLOPES
...                 BOX CARS and SLOPES
      §...                 PARABOLICAL
      e...                 (special: an exponential decay)

      t...                 TSF EDIT as a signal generator
(4)   +...                 LIN-COMB DATA


Select signals for a customized (hand-selected) predictions:

Bracket signal declarations with  BEGIN_HS ... END  commands:

BEGIN_HS [n1 [n2 ...]]
BEGIN_HS [n3 [n4 ...]]

BEGIN_HS [set-number [set-number...]]
BEGIN_HS [set-number [set-number...]]

See the example.

With this feature you can  generate a number of sets of predictions, set  n1 to set  nn , consisting of the selected effects
(tides, site effects or observed time series, multiplied by the estimated admittance).
In an upgrade of  urtapt.f , the maximum number of hand-selected sets is scalable
( include file  urtapuc.fp , parameter msxhs ).

Concerning CONDUMPS: Since we wouldn't like to output unused memory in the various dump operations, an option is available:
(type logical) in the namelist. If true, a slightly new dump control record is written and the output is limited to the necessary.
If false, the old convention, always 10 sets regardless of being used, is saved in the dump  )

The current maximum number of hand-selected sets is 20

What goes into the set is the segment of signal generating instructions that are bracketed by one pair of
BEGIN_HS ... BEGIN_HS  or a pair of  
With the former pair, one effect can be put into a number of sets.

If you don't specify a set-number, the system provides one by counting the number of  BEGIN_HS  statements.

However (!),
if qouts_i_sel (namelist) is set .true., the indicated constituents are excluded and all_other effects are included.

Detrending consists of a ramp signal and, optionally (if   q_allow_dc = .true. ), of a bias.
Since detrending may be selected automatically, you can decide where you include the trend
using the  opt_hand_select_trend  (namelist)  option: 

opt_hand_select_trend = 'string'
If the nk'th  character in the option string is  'T' or 'R' or '%', the ramp will be added to set nk.
If the nk'th  character in the option string is  'C' or 'B' or '%', the bias will be added to set nk.
If the option string starts with  'Own', a new set with only the trend is generated.
If the option string starts with  'No', no set will include the trend. Repeating:
opt_hand_select_trend = 'Own|No|string'
A partial residual can also be produced. using the  opt_hand_select_resid (namelist)  option:
If the  nk'th  column of the option string contains a 'T' a partial residual time series will be
produced (observation minus prediction).
The output time series are named
  rootname.phn1.ts  to  rootname.phnn.ts (predictions) 
to  rootname.rhnn.ts (residuals)

( rootname  is assigned by tsfile_names='rootname' in the namelist.)
                                              \ o *
Coding TIDE WAVE GROUPS:               ______*______

Identifier = (

(argument selector) = tide symbol

Argument selectors:

¨The set of astronomical arguments is
args = (n,tau,s,h,p,N,ps,pi/2)
The first selector specifies the spherical harmonic degree n. Therefore it must be a specific integer.
  For gravity tides, the smallest integer is 2, the largest allowed is determined by the tide potential. 
  In the case of solar radiation tides, the smallest integer is 1

The second selector specifies---implicitly---the spherical harmonic order (it multiplies the lunar orbital anomaly, tau). 
  Also here, only one specific integer should be given, since the response of any kind of instrument will be substantially
  different across the fundamental bands,
  0 for long-period, 1 for diurnal, 2 for semi-diurnal waves, etc.

For the eighth argument there is no selector, since the argument merely switches between sines, cosines and/or changes their signs.

The following applies to selectors 3 to 7:

 i   - an integer. Check with the tide potential development that the wave is represented at all.
 i:j - a specific range of integers, inclusive.
 i:  - from i to the largest argument multiplier present in the tide potential.
 :j  - from the smallest up to and including j
 *   - all waves

E.g.: (2,2,0,*) = M2      signifies the M2 band, all arg.numbers for arg h
      (2,1,1,-1,:0) = K1  the K1 band, all arg.numbers less/equal zero for arg p



Identifier = <

<argument selector) = tide symbol

This code must follow the wave group descriptor which is to be extended. Best described by example:

(2,2,2,-1:) = K2          designates K2 band
<2,2,3:)    = K2+         adds all tides beyond K2 to the K2 band


Identifier = (~

(~ frequency ) = symbol

Frequency in cyc/d

(~ 4.00d0 ) = S4 

Don't expect to easily resolve quarter-diurnal non-linear tides; there will be a set of terms with very similar frequencies, and since amplitudes and phases
won't be a straight affair, it's useless to combine them into a wavegroup. You simply have to calculate intermodulation
frequencies for candidates and add a number of  (~  -commands, with the risk to hit the roof of  the system matrix.

Coding BOXCARS,  SLOPES, PARABOLAS,  and EDECAY=Exponential Decay :

Identifiers =   =  % / P  §  e

Exponential decay: From a starting time given by a time-range parenthesis.
The amplitude will be estimated. (We cannot estimate a decay constant 'cause problem's nonlinear.)

e const [_/h_] ( {From date}{ To date} ) [ comment][ = NAME] = SYMBOL

const  - either a 1/e-decay time in hours or a decay parameter in 1/h

Box car, i.e. a time-limited offset that is to be estimated:

= ( {From date}{ To date} ) [ comment][ = NAME] = SYMBOL  

Simple slope, time-limited:

/ ( {From date}{ To date} ) [ comment][ = NAME] = SYMBOL

Simple parabola, time-limited:

P ( {From date}{ To date} ) [ comment][ = NAME] = SYMBOL

Box car and slope, time-limited:

% ( {From date}{ To date} ) [ comment][ = NAME] = SYMBOL

Parabolical with offset and slope, time-limited:

§ ( {From date}{ To date} ) [ comment][ = NAME] = SYMBOL

date   - 3 to 7 integers, space-delimited, specify
         year month day {hour minute second millsecond}

NAME   - a short string (4 chars) to identify the signal in the result protocol
         If not given, the first four characters of SYMBOL will be used

SYMBOL - a short string (typically 1 to 8 characters).
         Specify * as a wild-card; this event will apply to any
         data set.

From  and  To  intervals - In a sequence of commands, if the intervals are to be
         abutted, the dates should
not overlap. Leave a gap of less than the
         length of the time step!
         If the starting hour (t0) is exactly at a mid-time interval, note
         that the program rounds the time by subtracting one second.
         Example creating tightly connected box cars for a 24h series with t0=12.d0:

         % ( From 2011 02 23 00 00 00 000 To 2013 10 23 11 00 00 000 ) = %BS2
         % ( From 2013 10 24 00 00 00 000 To 2019 12 31 23 00 00 000 ) = %BS3

By default, the reference time of theses features except exponential is the centre time of the interval.
The default can be changed to the from-date with  namelist parameter q_tref_slopes_start=.true.
If NAME contains the character `/´ the reference time will be the center of the interval whatever the default is. 

NAME is significant for combining two and more BoxCars into a BoxTrain.
Specify the same name for disjunct boxcars, and one coefficient will be solved for the boxtrain.
Application: Absolute gravity campaigns where a setup is remeasured and the target quantity is the multi-project gravity offset.

NAMES should be unique to distinguish slopes and parabolas (and to allocate one array column for each term.
The commands % and § issue column generating subprocedures for two resp. three terms.
If   NAME  contains a  %-sign, it will be
replaced with `=´ for the boxcar, `/´ for the slope, and `J´ for the parabola.

If you want to create a boxtrain with slopes and you code identical names, as in
% ( From 2014 05 27 19 44 50 000 To 2014 05 29 04 53 25 000 ) = %SLO     "Onsala_AC_20140527a.drop.txt"
% ( From 2014 05 29 05 59 50 000 To 2014 05 30 05 53 25 000 ) = %SLO     "Onsala_AC_20140529a.drop.txt"
the following signals will be produced:
    One with a boxtrain named =SLO with the two boxcars; only one parameter will be solved for the train.
    Two slopes with the same name /SLO  for which two parameters will be solved (really??); this will make the result sheet slightly more difficult to read.
Yet, with  qpr_drift=.true. in the namelist,  TSF-EDIT commands will be produced in the result sheet with the appropriate time ranges.

Looking up Box cars and slopes in a SITE EVENTS table:
The site events file is a file consisting of Box car and Slopes definitions as introduced above.
The definition records are looked up using a   TARGET string and
comparing it to the SYMBOL:

= #file-unit (comment)  -> TARGET[ext]

file-unit - integer - file unit number, must have been opened in the file
                      opening block.

TARGET    - char*32 - is compared to site events file symbols.
                      All records with  target.eq.symbol  are taken.
                      Length of string in comparison is determined by the
                      shorter of the two.
                      Basically, the target would have four letters.

                      An extension would consist of two extra letters; for
                      uniqueness the first extra letter should be `_´
                      so that no ending of a four-letter-only symbol
                      would match a two-letter target extension.

                      Equivalence cases for TARGET's extension (last two characters):
                      _H  matches SYMBOL's last two characters _N _E _T _L
                          Use: events that affect horizontal components.
                      _V  matches  SYMBOL's last two characters _R _V _U
                          Use: ...vertical components.
                      _^  matches vertical and horizontal components
                          except those marked _|
                      _*  matches any component

= #22 (Site events) -> ALES

Note that you can use setenv to re-use the ins-file flexibly.
setenv SITE ONSA

= #22 (Site events) -> ${SITE}_H

There is more capability. Look into urtapu9.f  "If there is a # sign"  for details
Example from ~/tap/Upto_95/siteevents.dat:

= ( From 1994 09 11 To 1999 12 31 ) Box Car = AROT = ONSA_N
= ( From 1993 12 01 To 1993 12 31 ) Box Car = MADR
= ( From 1993 01 01 To 1994 01 31 ) Box Car = ASON = *

In this example, the event ending 1994 01 31 occurs in conjunction
with switching-on Anti-Spoofing on Feb 1, 1994 (we thus baptise it "ASON"). This event is likely
to affect all (TurboRogue) receivers. The wildcard `*´ accomplishes this.
If there are other receivers in the network unaffected of the anti-spoofing change,
you need to provide an explicit list for sensitive sites only, with the same event time - unfortunately.

The Antenna rotation at Onsala in September 1994 is marked with the symbol ONSA_N
If events are scanned with the target ${SITE}_H and $SITE happens to be ONSA ,
the event will be picked (since _H  matches both _N  and _E )

An event marked  ABCD_|  is a vertical-only event.
A target that is careful to not pick this event when a vertical component is solved should be coded



Identifier = t

t #iutse target = symbol

(You can code iutse as a string #nn where nn is a two-digit number.
You can specify it without the `#´ too.)
If you're going to exhaust the file unit range, use:

##openf-string for unit iutse
t #iutse target = symbol

The default signal is zero. The job will then fail at matrix inversion.
Spill space is granted (100%) and work space of 200% of the input data length.
At return, if the TSF EDIT operation has generated more data than the input data length, an error message is issued.
The situation is still benign as the overshoot will be ignored or overwritten with the next column filling.
If less data is generated, the trailing end will be filled with the DC-level of what has been generated.
A warning will be issued. This should be avoided. Try and fill the array to the end.


Designating OBS-DATA
= data input into signal matrix using Read_FuF through subroutine

Proc_f_on_rec ( ~/sas/p/urtapu13.f )

OBS! See Namelist parameter  q_del_miss_obs

Identifier = #

#iuninx, trg,rec_mrs,fmt,khms,itz,lddf [   E:tsf_editname] ]  = symbol [ -> EOF ] [ +M! ]

#iuninx, trg,rec_mrs,fmt,khms,itz,lddf [  E#:tsf_editname] ]  = symbol [ -> EOF ] [ +M! ]

#iuninx, trg,rec_mrs,fmt,khms,itz,lddf [   R:repairoption] ]  = symbol [ -> EOF ] [ +M! ]

#iuninx, trg,rec_mrs,fmt,khms,itz,[-9  AMP=value]             = symbol [ -> EOF ] [ +M! ]

Short version for BIN-files:

#iuninx [L:label ][ D:lddf ][   E:tsf_editname] ]  = symbol [ -> EOF ] 

iuninx       - integer  - file unit number

trg          - char*8   - target string, right-delimited by ]
                          'BIN' for binary data
                          '...] Rwd' to rewind file before read.

rec_mrs      - real*8   - missing record symbolic value. Significant only in the
                          case of ASCII-data. Binary data files keep an internal

fmt          - char**   - format string '(...)' for ASCII or 'U' for binary data or
                          'L:label]' in the case of 'BIN'ary multi-component
                          files. In that case this option must occur first.

khms         - integer  - expect 1,2,3 time fields (hours,minutes,seconds)
                          in the case of ASCII data.

itz          - integer  - time zone (0: UTC; 1: CET etc.)

lddf         - integer  - filter/scaling code value, c.f. below.

E:           - codeword - for calling TSF_EDIT (~/sas/p/tsfedit.f)
                          Commands will be read from Instruction file.
                          All tsf-edit sections must be collected in one file;
                          this concerns all file commands, i.e.
                          identifiers `#´ and `+´
E#:          - codeword - Commands will be read from the tsf_edit_file (Namelist)
tsfeditname  - char*16  - TSF_EDIT target string to locate the command section.

R:           - codeword - for calling REPAIR (~/sas/p/repair.f)
repairoption - char*4   - L to interpolate linear, else set to DCLevel

=            - codeword - identifier follows.
symbol       - char*4   - identifier "tag" for the data series used in results
                          table. Exactly one blank between = and tag.

-> EOF       - codeword - This time series participates in the Empirical Orthogonal
                          Function (EOF) setup. Makes it sensitive to the
                          q_auto_select_eu option.

+M!          - codeword - to allow missing records; they will be replaced with the DC-level
                          at the array contraction stage (urtapu2.f subr. CONMRS).

lddf         - integer  -  -9 < lddf < 0:  take data as it is; read no further line.

                          lddf = -9:  special value, see below
                          lddf =  0:  next line is (READ (4,*)) ...
                                      sf [PCAL]      where sf = ScaleFactor
                          lddf >  0:  next lines are (READ (4,*)) ...
                                      sf, nfm,nfp,symm [PCAL]
                              sf          - real*8  - scale factor for this data set.
                              PCAL        - codeword- use scale factor as conversion of units.
                                                      [input file]/[pressure units].
                              nfm,nfp     - integer - index range of the filter that follows;  
                              symm        - char*1  - filter symmetry code ('A' | ' ' | 'S')  
                              lddf=-9:                special value
                                                      AMP=value   or  AMP=name
                   must be specified later on the same line
                   This is for multiplying the time-series with a factor.

value       - real*8  - a numerical value: scale factor

        - char    - a reserved name, use  fgrep -i recamp ~/tap/m/urtap.f ~/sas/p/*.f

                        to explore what is available.
                        There are PTX and PTY associated with namelist variables ptfx and ptfy

#23,'BIN]',-9999.9d0,' ',0,0,-1 = AplO
for air pressure loading, reading binary data. (This is almost the simplest possible example, though there are redundencies.)
We take the data as is; khms is ignored in binary data, and itz is unchanged. The last `-1´ (lddf) means we don't filter nor rescale. The -9999.9d0 as the missing-record symbol does not apply in binary data sets.

The usual way to refer to a file is by including it in the open file block. However, if many files need to be loaded into the observation array, this method will exhaust the set of possible unit numbers. There is the following way out:


i.e. add a ##-command line before the #-command. The string following the ## is passed to the openf procedure (~/util/p/openf.f) via the  entry openfs. After proc_f_on_rec (~/sas/p/urtapu13.f) the file will be closed. Example:

Tide Bands
(1:2,0,0,-1:1)  = Sa
(1:2,0,0,2,*)   = Ssa
<1:2,0,0,-2,*)  = Ssar
(1:2,0,0,3,*)   = Sta
<1:2,0,0,-3,*)  = Star
(1:2,0,0,4,*)   = Sqa
<1:2,0,0,-4,*)  = Sqar
##23 O ~/Oload/smhi/apl/SUND_apl.tsf
# 23,'TOP]',-9999.9d0,'(t6,a9,1x,i2,t18,f10.3)',1,0,-1,-1.0,24.0  = aprL
= #22 (Site events file, look for) -> SUND_U
End Tide Bands
This example reads Sundsvall air pressure loading time series (6h samp) from unit23 into the signal model array (24h samp). The string passed to openf is "23 O ~/Oload/smhi/apl/SUND_apl.tsf"

Side-effect: With this method of file opening the file name does not become available through the file name common-area (/COPENF/),  yet some of the subfunctions might like to use it for commentary purpose (displays, protocol). An "explain"-feature is iunder development (June 2015); anything enclosed in " " will be saved as an explanatory comment to the column / solution coefficient.

## is synonymous with #_   #:   #<   #.   #{    #(

Here's how to input a series from a multichannel file:

##23 # /geo/hgs/geoc/jpl/
# 23 L:X_VAL D:0,0,0 R:.] = GEOX
##23 # /geo/hgs/geoc/jpl/
# 23 L:Y_VAL D:0,0,0 R:.] = GEOY

##23 # /geo/hgs/geoc/jpl/
# 23 L:Z_VAL D:0,0,0 R:.] = GEOZ


We access the data by labels X_VAL etc. (there are also X_SIG etc.). We rescale JPL-geocenter motion from millimeters to meters, and repair missing data by substituting DC-level (R:.]). The file is closed automatically after each read, so we need to open it everytime. We could have chosen a unique file unit and opened the file from the leading file-opening block; but then
we need to add a rewind option ( Rwd) in the format string: 'L:X_VAL] Rwd' etc.

Outlier detection and outlier editing
If editing of outliers is desired, the job must be run two or more times. In the first run the outliers are detected, then applied, new outliers are detected and added, etc.

To edit outliers, one of two  TSF EDIT options is used. A suitable block would be

DEL From 1993 11 15 0 0 0 0 To 1994 01 10 0 0 0 0
DEL From 1995 02 04 0 0 0 0 To 1995 02 04 23 0 0 0
CONT 31 O d/LEKS.ra.tse

For early editing, the Namelist &param would contain
 tsf_edit_name='LEKS]', q_tsf_edit_y=.true.

For backward compatibility the parameter  q_tsf_edit is still recognized as the early-editing alternative.

For late editing,
 tsf_late_edit='LEKS]', q_tsf_edit_late=.true.      [sic!]

See remarks on early and late editing below in this subsection.

Both reading and writing outlier records is possible in the same job.
Thus, re-running the procedure builds up a growing outlier file.
Iteration can be performed until the protocol (unit 6) reads
 <Main-->>> No outliers detected.

File unit numbers:
The  CONT  input unit number 31 in the  TSF EDIT  section is used merely because we know the unit is available.

The  iun_outliers  unit number 31 implies a special function though: The file name is generated from the  tsfile_names  parameter (extension .tse  is added), and the file is opened for append. In the case of a unit number different from 31 the unit must be connected to a file declared in the top open-file block, and overwrite or append can be controlled with the options field.

Outlier criterion:
A positive number detects an outlier as a value differing from the mean by the given absolute amount.
A negative number detects an outlier as a value differing from the mean by the specified factor of standard deviation.

To write outliers, the Namelist &param would contain
 tsfile_names='d/LEKS.ra '
 iun_outliers=31, outlier_criterion=-5.0

Early and late editing:
"Early" means after reading, scaling, and trimming the observations and before computing the tide potential.
"Late" means after filtering and right before gap contraction.
Filtering before editing will retain more samples; however, these could eventually be compromised due to
far-lag filter responses. On the other hand, long filters can devour a lot of data if the outlier editing results
in many and tightly placed gaps.

OBS! Late editing implies that the residuals are already output before the TSF EDIT instructions are applied.
They would typically contain  DEL  commands for outliers. If the cleaned residual is to be plotted, e.g. using tsd,
then use the urtap-instruction file
   tsd -n #n... -E urtap-xxx.ins CLEAN file.wr.ts
(if  CLEAN  is the locator for that section in that ins-file).

There is another advantage with late editing: The residual can be used subsequently with our without the
outliers. Use
    tslist file.ra.ts -Eurtap-xxx.tse,CLEAN -I -o file-ed.ra.ts

There is more information in the section dealing with the Namelist.


Customized residual.

(This is a pre-runner of the "hand selected residual". Since it's more easy to read, I have left the paragraph here to stay.)
The system computes tidal and nontidal residuals, and total residual. This leaves room for a wish. To request a tailored residual the user can surround the wave group instructions with BEGIN - END blocks. Thus,

Tide Bands
(1:2,0,0,-1:1)  = Sa
(1:2,0,0,2,*)   = Ssa
<1:2,0,0,-2,*)  = Ssar
(1:2,0,0,3,*)   = Sta
<1:2,0,0,-3,*)  = Star
(1:2,0,0,4,*)   = Sqa
<1:2,0,0,-4,*)  = Sqar
= #22 (Site events file, look for) -> VANE_U
End Tide Bands
devises a partial residual which will be composed of
r = y - A x'

where y are the observations and x' are the coefficients of the signals in the  B E -block. In the example it means that the jumps from the site events file are removed from the observations.

More than one such block can be specified.

If the feature is in use, the program will output the following files - time series of  A x'
tsfile_name.rh.ts - time series of  r

In order to do the inverse of the above, i.e. exclude the indicated components, specify in the namelist


Designating LIN-COMB DATA, reduction of data to be fitted,
by linear combination of data which is input using  Read_FuF.

Identifier = +

Short version for BIN-files:

+iuninx [L:label ][ D:lddf[,vadd] ][ E:tsf_editname] ]

Long version for ASCII-files:
+iuninx, trg, rec_mrs, fmt, khms, itz[, lddf[, vadd ]] [ E:tsf_editname]  ]  

Again, lddf declares whether scale factor or filter information follows.

A simple scale factor is applied using  lddf = -1  and  vadd = scale
Default lddf,vadd  is -1,-1.0

vadd   - real*8 - +1.0d0 | -1.0d0
                  add or subtract the data

                   to or from     the observations

Other parameters c.f. OBS-DATA BANDS


Y <- Y + vadd × ( Filter * R )

The filter may have just one coefficient, ddf(0), and lddf=-9 signals that
AMP=#v  or  SCA=#v
(appended to the instruction string) sets its value. SCA will set dff(0)=1/v
In that case, note the multiplication of vadd and dff(0) 

The  AMP / SCA  versions of scaling can be combined with environment parameters.

Note that the data must be in synch with the rest. If the start cannot be adjusted, urtapt will stop.
If the series is too short, there are three options to continue processing: cut, delete, prolongate.
See namelist parameter opt_missalign

Mixed tides can be analyzed (CAUSE='TGP|SRT'). Two batches of tide band
information must be given (sections (1),(2) (and (3)) are doubled,
sections (4) and (5) append the second batch.

        Tide Bands
        * Switch to SRT
        (1:2,1,1,*) = S1
        (2,2,2,*) = S2
        End Tide Bands

Input from files:

IUIN      - File unit for input
IUMC      - For multi-component (multi-channel) data:
            Run first job with first component input from binary multi-component
            obs data file IUIN. Successive jobs: Process dumped data (QGet_Dump
            =.t.) and read data from unit IUMC. User is prompted to select the
            IUMC+1 is used for temporary storage.

41 (42)   - Prediction error filter (PEF) bank. Use unit 42 for read-only.
            Open both to read from 42 and write to 41.
            Usage: 42 for an a-priori filter, 41 to re-iterate it.

82        - Previous solution parameters, *.evs file (out unit 12)
            for iterating solution parameters, especially in the case of
            a priori model covariance data.

            The following files are used for incorporation of e.g. air pressure
            corrections in the tide potential. Processed under CAUSE='TGP+SRT'

91        - An SRT-job's *.dmp file (out unit 11 under CAUSE='SRT')
92        - An SRT-job's *.trs file ( "   "   13   "     "      " )

Input data can be read from a multi-column file on unit 21. For this purpose, namelist input should include
 TRG='BIN', FMT='L:label '

Since environment parameters are resolved in FMT outside a parenthesis, one can specify
 TRG='BIN', FMT='L:${site}|${tapcomp} '
to retrieve a label 'ONSA_________RAD' after
  setenv site ONSA
  setenv tapcomp RAD

Output to files:

File unit reserved for this and that: 9

Files that must be presented in the open-file block at the top of the ins file:

IUN_RESULTS - RSLIST result sheet:                                        *.prl-file
              IUN_RESULTS is a namelist parameter [default=14].

11          - dump of observations, signal model, tide potential etc.     *.dmp-file
              Reprocessing (change number of eigenvalues, transm-mode,
              length of noise whitening filter, iterate further):
              Set  QGet_Dump=.true.

Optional files the open status of which is inquired; if open, write ...

41          - Prediction Error Filter bank of last iteration              *.pef-file BIN
12          - eigenvalues, -vectors and solution parameters:              *.evs-file BIN
13          - tide band information and response parameters:              *.trs-file ASCII
15          - Variance/Covariance matrix                                  *.vcv-file BIN
15            alternative with  qout_vcv_bin=.false.         
             *.vcv-file ASCII   
19          - Solution array and design matrix                            *.gmx-file BIN

       - Abridged results of multi-component analysis output.

              IUMCR is a namelist parameter [default=16].

IUN_SBSEL_OUT - Selected columns of the signal array after filtering but
              before condensing (conmrs), output to labeled MC-file       *.mc        BIN                
IUN_CONDUMP - Input of merged dumps (from urtapm) or output of dumps
              for merging (to urtapm).
IUN_CONDUMP is a namelist
              parameter, and so is OPT_CONDUMP specifying
              I/O-instructions. [default IUN_CONDUMP=0, no I/O]           *.cdmp-file BIN

18          - With condump input, the header lines are written to         *.lst-file  ASCII
              this file.               

Eigenvectors in data space:
- controlled with namelist parameters ndump_eu (the amount of
              vectors belonging to the greatest eigenvalues)
fndump_eu (file name). 

Time series output (predictions, residuals) is via file unit 31.

File names are compiled by the program, using namelist parameter   TSFILE_NAMES
and adding the following extensions:
.tp.ts          "tide potential" = predicted time series using unity tide coefficients
                (if qouts_tp=.true.)
.pn.ts          predicted non-tidal (box cars, slopes, observed ancillary)
.rn.ts          residual  non-tidal
.pt.ts          predicted tidal     (band command codes "(" and "<" )
.rt.ts          tidal residual
.ph##.ts        predicted custom signal (if declared in band section between BEGIN and END)
.rh##.ts        custom residual. The ## is replaced with the Hand-selection set number.
.pa.ts          predicted "all", pn + pt
.ra.ts          residual w.r.t.  pn + pt  (with namelist q_ra_with_outliers=.true., with outliers)
.fa.ts          filtered residual
.sr.ts          smoothed prediction (currently with a deterring problem when using a diff-filter)
.dw.ts          data weights (if qouts_weights=.true.)
.ds.ts          sigmas = reciprocal weights (if qouts_sigmas=.true.)
.wr.ts          filtered and weighted residual = the actual residual at the
                least-squares solution stage.
.jd.ts          in conjunction with CONDUMP, the modified Julian dates.

.po.ts          integer arrays; positions of outliers in conjunction with CONDUMP.
                This file is re-input in successive runs, and new records will be appended
                if additional outliers are detected.

.tse            ASCII: outlier records if IUN_OUTLIERS=31

The .ts-files are binary

File name compilation involves replacement of environment parameters. E.g. the following can be specified in the namelist:
 TSFILE_NAMES = '${site}${tapcomp} '

Environment parameters 

URTAP_LDATE  - if namelist &param specifies LDATE(1)=-1, $URTAP_LDATE will be used.

URTAP_NEXTRA - if namelist &param specifies NX_EXTRA=-1, $URTAP_NEXTRA will be used.   

URTAP_NDATA  - length of data array for analysis. If not given, Namelist parameter
               NData_Use is taken.
               You may use
               setenv URTAP_NDATA `chck-tsends urtap.ins`

Note: Truncation of long input series must allocate space for reading; set NX_EXTRA to a large enough value
(at least the length of the longest overshoot). Unfortunately, urtap cannot give a hint.

In a few (currently four) namelist variables (of type character) you can introduce environment parameters and
set values before the program is called. The reference syntax is different from csh, 

E.g.    tsfile_names = 'o/g${URTAPFD}-${URTAPDT:[1h]} '


Source code PARAMETERS

Array size

       NYD    - max. time series length.

       NBANDS - number of signal components
                (2 for every tidal wave group).

       MXSPEF - maximum order of PEFilters.



Parameter               Meaning                                      [ default ]

_____________________ EXEC CONTROL _____________________________________________

QBATCH                - logical - Avoid all prompts and graphics.

                      If QBATCH=.T.
                      the following specifications are substituted at
                      the corresponding prompts:

use_window_residual   - char*1  - "Change status window on the residual"
                                  recommended:  'F'  (ascertain false)     ['N']
                                  Status at wake-up: Q_window_residual

use_downw_residual    - char*1  - legal values: 'Y' and 'N'
                                  to downweight the residual without
                                  a window; overrides
use_window_residual  ['?']
                                  (Was introduced since use of windows is

use_eigenvalues       - char*16 - User response to "Select No.of eigenvalues"
                                  prompter in interactive mode

                                   e.g.'W 6.0'
                                  'G' for Go! is automatically added.
                                  Possible responses:
                                  W #w   - cut away the smaller Ev's
                                           beyond 10-w of the greatest.
                                  D #n   - delete the n smallest.
                                  I #n   - include the n greatest.
                                  A      - include all.                    ['A']

use_solution_options  - char*16 - "Solution options"
                                  Possible entries:  +U +L +T +S +E R
                                  +E     - EOF-analysis
                                  +U  +L - linear alt. Log scale (in plot)
                                  +T     - Transform with model var/covar
                     Have to check this!
                     In urtap/urtapt lines 887ff +U seems to have two meanings. Check Subr. Set_Solve_Mode

                                  or blank                                 [' ']

Q_Get_Dump  - logical - Rerun previous job; take data from dump (unit 11)    [F]
N_Old_Dump  - integer - Dump format of older versions (-1 -2 ...)
                          for compatibility                                  [0]
_____________________ DUMPS ____________________________________________________
Here, a new kind of processing is introduced. Keyword CONDUMP.
Dump the condensed Design matrix and the associated arrays such that a suite of dumps can be concatenated,
and solutions of campaigns can be produced where campaigns might be widely separated in time.
A series of dumps can be merged with program urtapm . Application is primarily intended for multi-campaign
analysis of Absolute gravity measurements.

Repeated execution of urtapt will iterate outlier editing using a mechanism that differs from the tsfedit method1):
We poke a small value into the design matrix and the observation array.  
This may have undesired side effects, which can be avoided by setting the criterion to a very large value;
or by deleting the outlier file produced: $tsfile_names//'po.ts'
The outlier detection is still done in the standard way (outlier_criterion, q_outliers_filtered)

iun_condump - integer - if > 0 dump to this file unit (must be opened in
                        open file-block

opt_condump - char    - 'GET' or 'GET UNW' - import and, with 'UNW' undo
                        weighting. UNW is discouraged; urtapm ought to
                        produce reasonable weights, in the first case
                        those from the original drop files. 
                        'PUT' export the design matrix and stuff.
                        String may include '+DBG', a debug option for
                        get or put /CRECB/ (esp. concerning hand_select)

            - logical - condense (and output) the weights.
                        set .true. if a dump file unit > 1
                        is declared and q_downw is .true.              [.false.]

rwpo        - real*8  - In outlier editing, the small value that will
                        replace elements in design matrix and
                        observation array at outlier positions.        [1.0d-10]

1) the reason being that the rows of the design matrix cannot be deleted without making the feature table inconsistent.
A tse-file will still be produced. This file can be concatenated with single-campaign outlier files in re-processing of the
whole multi-campaign solution from the beginning (urtapt on the affected campaigns, re-run of urtapm, new feature table, updated urtapm-tse files ... )

            - integer - File unit for dump of selected columns of the signal array
                        right after filtering and before condensation       [-1]
            - integer - The band numbers (not the column numbers, sorry!)
                        Obtain them from the log file or calculate:
                        Every harmonic band will be output in cos and sin;
                        The prl-file shows the odd numbers n <= 2*NTB-1 of
                        the columns for them. So the first column that a
                        non-harmonic signal is assigned to is 2*NTB+1.  [-1,...]
           - char*128 - Instead of coding the numbers into
k_sbsel_out you can
                        write expressions like 'j1-j2,j3-j4,j5,j6' into
                        this string                                        [' ']
                        This option, if non-blank, overrides any
                        assignment of

ulab       - char*16  - (only the first 10 are significant)
                        The output MC-labels' first characters,       ['GARRAY']
                        the remainder is filled with the band symbol.

_____________________ TIDES and other signals; DC ______________________________

THD_file     - char*64 - Tide potential file     ['/home/hgs/tap/d/etcpot.dat ']

CAUSE        - char*8  - 'TGP'     - tide generating potential
                       - 'SRT'     - solar radiation tide
                       - 'TGP+SRT' - air pressure perturbed tides
                       - 'TGP|SRT' - TGP and SRT side-by-side         
                       - not 'TGP'
                         nor 'SRT' - no tide will be generated;
                                     observations may be specified in
                                     the `Tide Bands´ block.             ['TGP']
             - char*32 - Darwin symbol or argument numbers for a
                         demodulating frequency.
                         E.g. '2,1,1,-1,0,0,0,0'                           [' ']
                         (under CAUSE='TGP' and/or 'SRT')

             - real*8  - Assume that the effective tide in the
                         signal model is the average tide during
                         the specified interval (in hours)               [0.0d0]
                         (under CAUSE='TGP' and/or 'SRT')

Iun_FCP      - integer - >0: Apply a filter on tide generating potential; 
                         read from file unit iun_fcp                         [0]
Q_allow_DC   - logical - .true.: The setup of a DC tide band to solve for
                         a bias is allowed (the weighting and filtering
                         procedures might leave a nonzero mean). Urtap will
                         have a look and decide.                             [T]

Q_always_DC  - logical - .true.: Always analyse a DC signal.                 [F]

Q_weights_DC - logical - .true.: The DC-level in the signal model will be
                         subjected to weighting like all other channels      [F]

Q_copy_DC    - logical - .true.: Take the DC-level of the input signal and
                         copy it into the predicted time series              [F]

DC_Level_Input  real*8 - Add this value to the detected DC-level of
                         the input signal                                 [0.d0]


             - logical - .true.:  Slope signals start with zero at
                                  the "From" -time and end with 1.0
                         .false.: Slope signals run from -1. to +1.          [F]

Q_RemDC_WMean¨- logical - Remove a DC-level from the input data using a
                          weighted mean                                      [F]

Q_FilterB_RemDC, Q_Conmrs_RemDC
              - logical - Remove DC-level in design matrix by REMDC
                          (not weighting) in
                          1. subroutine filter_bands
                          2. subroutine conmrs      

_____________________ SITE AND TIME ____________________________________________
XSITE,YSITE  - real*8  - Site coordinates, geographic, E-long., N-lat.       [-]
Idate(3)     - integer - Requested date of epoch for reference. A negative
                         Idate(1) uses the first date of the input file.
                         Used with Q_Cut_At_Epoch=.true. will truncate
                         the input time series.                           [3*-9]

Start_h      - real*8  - Requested hour after epoch for first sample      [0.d0]
                         Used with Q_Exact_Date the reading routine will
                         be called with t0=start_h. In that case, DT
                         must also be specified.
                         Used with
Q_Cut_At_Epoch=.true. will truncate
                         the input time series. 

Ldate(3)     - integer - Requested ending date,
                         Ldate(1)<1900: as much as available              [3*-9]
                         Ldate(1)=-1: Read from environment URTAP_LDATE

DT           - real*8  - Enforced sampling interval (e.g. to undersample)
                         Zero value: DT is determined from input data,
                         to be preferred with BIN input.
          If DECIMATE is used on input data (discouraged!), a nonzero
value leads to misrepresentation in results                     [0.d0] 
                         (??? Parameter setting for the result sheet has
                         taken place before tsfedit, even before the
                         "early" call ??? - should be fixed)


DTUT         - real*8  - D(TDT-UT) in [s]; IERS-Bull.B: D(TAI-UT) + 32.7  
                         Is obsolete unless SRT is used                   [57.0]
IALTG        - integer - Global reference for grav. and site radius.         [2]

Q_Cut_at_Epoch, Q_Exact_Date                                  [.false., .false.]

UTC_file     - char*64 - File name for leap seconds

Q_Inject_Leaps - log   - Inject leap seconds in tide synthesis          [.true.]

_____________________ INPUT DATA ______________________________________________________
IUin         - integer - log.unit for data to be analyzed, Read_Fm_D        [21]
IUmc         - integer - multi-component data file under QGet_Dump=.true.
                         A value less than 1: use the data dumped on file
                         unit 11                                            [-1]
TRG,...,ITZ            - parameters for Read_Fuf. Notice that...
FMT          - char*64 - may include environment parameters after the
                         last closing parenthesis.
L_Data       - integer - >0: Max.duration for reading procedure.              
                         <0: Truncate by this much
                         =0: Use all                                         [0]

Ndata_Use    - integer - Use input data up to this length,
                         =0: use all                                         [0]

NX_Extra     - integer - Additional space for reading #-bands             [1000]
                         =-1: read from environment URTAP_NEXTRA

Data_Mrs     - real*8  - Symb.value for missing data used internally   [-9999.9]
TST          - real*8  - Test accuracy for Data_Mrs / Rec_Mrs              [1.0]
Scale_Y      - real*8  - Scale factor, multiply with input data            [1.0]
Y_Units      - char*16 - units for data under Display_Ts                ['[m]$']
                       Editing input data:

Q_TSF_Edit_Y     - log - Call TSF_EDIT to work on the input time series      [F]
TSF_Edit_Name -char*32 - name tag of the "early" TSF EDIT block               []
TSF_Edit_File -char*64 - Use this file with this ...
IUN_TSF_Edit - integer - unit number                      [-1: instruction file]

D1_Rec      - char*128 - Supplements first input record to set a date.
                         Only the date part is significant.
                         Used in a call to Read_Fm_D1.                  ['NONE']

             - logical - Missing observations in the design matrix can
                         either be treated as no signal (zero) or
                         by removal of the affected rows of the design
                         matrix ("contraction").
                         A value .true. will ask for the latter        [.false.]
                         Note that missing obervations in the
                         r.h.s.-series will always lead to contraction.
                         Use tsfedit-options if you e.g. want to interpolate
                         missing data.

             - char*1  - How to deal with lin-comb-data that is too short.
                         If less than 10% is missing at the end, there are
                         three options:
                         'C' - cut; 'D' - delete; 'P' - prolongate
                         'Cut' adjusts all series' lengths to the minimum;
                         else the series is extended:
                         'delete' injects missing-record symbols,
                         'prolongate' repeats the last value.
                         Any other code: Processing will stop.             ['S']    

_____________________ DATA: NOISE (PRE-)WHITENING ______________________________
DFF(0:L_Dff) - real*8  - Difference filter, L_Dff = 0 | 1, used for
                         PRE-whitening (OBS! the default)
                         If L_Dff=1, this filter is always involved
                         (LPef and filter F() are alternately staged
                         after DFF)                                   [1.0,-1.0]

LPef         - integer - Use a filter (from file unit 42) of this order
                         in first program cycle. Overrides specification of
                         F() (c.f. below).
                         LPEF <= 0: use only the DFF-filter                 [-1]

F(0:L_Filter) -real*8  - Filter, must be nontrivial if L_Dff=0             [1.0]
                         Will be refined using Burg-proc.
S_Filter     - char*1  - 'S' - symmetric, 'A' - asymmetric                 ['A']

QBurg        - logical - Produce a new PEF                                   [F]
MXLPef       - integer - Maximum order for Burg PEF determination.          [20]      

             - logical - Find the longest uninterrupted section
                         as the input to the MEM analysis                    [T]
i1_Burg      - integer - First sample for Burg. If the first section of
                         valid data is too short, Burg will cause program
                         to quit. In the case of fragmented data check
                         with time series display of residual to find a
                         good starting position.                             [0]
m_Burg       - integer - If m_Burg > 0 the time series length from which
                         the PEF is determined will be limited to m_Burg    [-1]
                         m_Burg < 0 : take all available data.
Ls_MEM      - integer - Length of MEM spectrum                       [Ls_Psp]
                         N.B. Under periodogram, ls_mem should be given
                         a reasonable value (smaller than Ls_Psp)
Quiet_Burg   - logical - Shut up, J.P.!                                      [F]
Detrend      - char*1  - Detrend the data ('Y'|'N'|' ')
                         'Y': Force detrending.

                         ' ': Automatic feature: Detrend if DFF has
                              nonzero response at freq=0

                         'N': do not detrend even though DFF...            [' ']  

Q_Downw      - logical - Process a Downweight block in *.ins file            [F]

             - char*64 - (under developement)
Controls the removal
                         of DC-levels in the design matrix after weighting.
                         Presently only 'AUTO' is recognized.
                         Without AUTO, control is possible with the
                         options in the Downweight block, especially
                         the command lines
                         SELECT KEEP (or SELECT REMOVE).
                         'AUTO' - DC-levels are removed from tides and
                                  data columns input from files         ['AUTO']

NRJD         - integer - Reject diurnals & overtones (S1,S2,...)            [-1]

_________________________ DC-LEVELS and SLOPES _________________________________

Q_RemDC_WMean¨- logical - Remove a DC-level from the input data using a
                          weighted mean                                      [F]

Q_FilterB_RemDC, Q_Conmrs_RemDC
              - logical - Remove DC-level in design matrix by REMDC
                          (not weighting) in
                          1. subroutine filter_bands
                          2. subroutine conmrs      
              - logical - Analyse the signal matrix's columns for slopes
                          and report.                                        [F] 

Tsys_Slopes   - char*64 - Exempt symbols with those characters, `.´ as
                          a wildcard. E.g. '.BS.' for a drift model that
                          is instructed with %BS1 =BS2 etc. Use `,´ as a
                          delimiter, e.g. .BS.,.DEC to exclude box cars,
                          slopes, and exponential decays.               ['....']

______________________ OUTLIER EDITING _______________________________________
                                       Producing DEL commands_________________
- integer - file unit to which outlier records are written.
             - real*8    < -1: outliers are values |y| > |c|*WRMS
                               where WRMS is the weighted residual RMS.
                         = 0:  ignore outliers.
                         > 0:  outliers are values |y| > c                 [0]

             - logical -  .true.: Scan the filtered and weighted residual
                         .false.: Scan the raw residual                    [T]

             - logical -  .true.: Print details how the DEL instructions
                                  group outliers together                  [F]
                                        Input and applying DEL commands_______

Q_TSF_Edit   - logical - .true. required to consider scanning the          [F]
                         ins-file for a TSF EDIT block with the target
                         string declared by ...
TSF_Edit_Name -char*32 - the block must then contain a command
                         CONT 31 O the-file-with-the-del-commands
                         If this file is identical to tsfile_names + .tse
tsfile_names is the value of the very
                         namelist parameter, then a cycle urtap runs
                         will accumulate DEL commands until no more
                         outliers are found.    

             - logical - Call TSF_EDIT to work on the time series before
                         it is scanned for outliers                        [F]

TSF_Late_Edit -char*32 - name tag of the TSF EDIT block                     []
                         The DEL instructions are expected in the same file
                         as the "early" editing instructions

_____________________ EIGENVALUE ANALYSIS ________________________________________

           - logical - .true.:  Demands all eigenvalues > 0
                        (the computed ones are squared; however,
                        numerical precision might cause negative
                        values). If |ev(1)/ev(j)| > 102w , the sign
                        will be reversed, else the program will stop.
                        (the spectral width w can be specified
                        at the prompter or, in batch mode, through
                        namelist parameter use_eigenvalues).

                        .false.: Negative ev's will be sign-reversed.

q_unimag    - logical - Make signal array columns unity magnitude. This will
                        result in a more balanced eigenvalue spectrum.      [F]
                        Obs! Consequences on eigenvalue dump and error-
                             plotting pgm grdcvem not tested yet !

eof_grades() - real*8 - grades for explicit eigenvalue sorting. Is used when
                        q_user_select_eu is .false. Grades are to assess the
                        amount of information that an eigenvector contributes.

                         Association of eof_grades with eval is in reverse
                        order, i.e. eof_grades(1) is for the smallest eigen-
                        vector etc.
                        Remember: The eigenvalues are already sorted by
                        decreasing magnitude.

                        Assigning a low grade, e.g. eof_grades(j)=j will push
                        the (n-j+1)'th eigenvalue to the end of the spectrum;
                        there it can conveniently be deleted at the
                        "select eigenvalues"-prompter. By default, a value of
                        zero will assign  eof_grades(n-j+1)=2*n+1-j. If no
                        eof_grades is assigned to deviate from this rule
                        there will be no change in the ordered set.
                        The use of this option requires knowlegde of the
                        Dimension eof_grades = nbands [50]      [nbands*0.0d0]

            - logical - Interactive prompting for grades

                        specification, showing the data-space eigenvectors
                        in graphic display. Is an alternative to  uval()          uval?! maybe info in urtapu22.f
                        especially when the user does not yet know the            Or maybe we mean the eof_grades?
                        eigenvalues/-vectors                              [F]

            - logical - under qbatch=.true., this option will utilize eigen-

                        vector projections to find out the deletable EOF's
                        automatically. The time-series that are the basis for
                        EOF must be marked with ' -> EOF ', see section on
                        time-series input and EOF guidelines.
            - real*8  - under q_auto_select_eu=.true., grace margin to admit

                        eigenvalues into EOF space. Normally, the criterion
                        is set at the value 0.5. The program iterates, but it
                        needs a starting value that returns rather too many
                        eigenvalues than too few.                        [0.0]
                        Try auto_select_eu_margin=0.1

_____________________ OUTPUT _________________________________________________
TSFile_names - char*64 - Output predictions and residuals va file unit 31
                         using  TSFile_names  as a name root to which
                         suffixes will be added.
                         '?' or illegal characters - don't output.       ['?']
IUN_results  - integer - Subr. RSLIST result sheet                        [14]

             - logical - The residual and the filtered residual will
                         contain the outliers (i.e. the observations
                         before they are processed by tsf_edit)

Q_ra_with_DC - logical - The residual .ra. will receive the DC-level
that was removed due to q_remdc_wmean=.true. [.true.]

Q_res_remdc  - logical - All residuals are produced DC-free          [.false.]

opt_hand_select_trend  - char*20                                ['..........']
_hand_select_resid  - char*20                                ['..........']
qouts_i_sel  - logical - "inverse selection"                         [.false.]
q_hs_as_needed           - logical -                                 [.false.]
                         ... see hand selected predictions

qout_vcv_corr  logical - Output on unit 15 if open:
                         .false. : Variance-covariance matrix
                         .true.  : Correlation                      
qout_vcv_bin - logical   .false. : ASCII output else BIN

ndump_eu     - integer - The amount of eigenvectors in data space belonging
                         to the greatest eigenvalues for output to ...    [-1]
fndump_eu    - char*80 - ... MC-file with this name           ['']

_____________________ RESULT SHEET____________________________________________
Q_Hindsight  - logical - If no uncertainty information is available,
                         the uncertainty of the estimated parameters
                         are inferred using the residual RMS          [.true.]

Q_Gravity    - logical - Tidal gravity. Affects rslist results only.
                         (will show delta-factors)                   [.false.]

Delta_system  - char   - WAHR  or  DEHANT (the latter is the Melchior-
                         Ducarme-Bruxelles preference to divide the observed
                         by the local astronomical gravity. We do that 
                         in a different way, see Delta system         ['WAHR']
p_units       - char   - Either '[m]$' (=default) or '[nm/s^2]$'
                         in the latter case, show the rigid-earth
                         gravity tide in the Amplitude column (else
                         potential elevation in [m]).                 ['[m]$']

Q_Inertial_Acc  - log  - In rslist results, in the delta-factors, the
                         effect of inertial acceleratiom is removed
                         if .true.                                   [.false.]
_Inertial_Acc - char - The file name from which the inertial terms
                         are read. They are produced with e.g.
                         deltafactm -iac -N20 -A -O d/Onsala-delta-inertial.dat
                                                                         [' ']
                    OBS! Standard is to include inertial terms in the
                         delta-factors, i.e. to leave the estimated deltas
                         unreduced (our default). Use deltafactm to produce
                         modelled deltas and ~/Ttide/SCG/tide-phasors
                         to make plots with ocean loading phasors etc. 

IALTG        - integer - Gravity value and radius at the site -
                         computation option:
                            *0: Gravity, Equatorial reference;
                            *1: Gravity, IERS-definition;
                            *2: Gravity, Local reference.
                         0* 1*: Radius,  Equatorial reference
                            2*: Radius,  Local reference                 [2]      NEW IN LSPOT.F !

GE_IMPOSE    - real*8  - Value to supersede the computed gravity value;
RE_IMPOSE    - real*8  - Value to supersede the computed radius value
                         (both computations take place in LSPOT.F; the
                         _IMPOSE parameters, if given here, will be used).
                         IALTG will not take action if reasonable
                         _IMPOSE values are given                  [0.d0,0.d0]   

Q_Ocean      - logical - set oceanographic phase convention                [F]
IUN_Fcp      - integer - File unit number, filter coefficients for

                            tide potential (calling F_On_Pot )            [-1]
IUmcr        - integer - MC results output unit                           [16]
Sp_Zero      - real*8  - Spectrum cut level (suppress below)           [1.d-8]
List_TGP     - integer - Amount of major tides to print                   [-1]

             - char*8  - Representation of slopes in result sheet,

                         e.g. units per year.
                         Possible values: 'year' 'day' 'cty'
                         cty = centuries                              ['year']

fmt_results  - char*64 - Format directives, (urtapu5.f variable),           examples
                         A:fmt  Numeric values (FNUMAMP)              A:1p,d12.4   A:0p,f12.6
                         D:str  Units degrees  (FDEG)                 D:5h[deg]    D:3h[o]    
                         S+ D:3h[^]  Special characters for degrees
                                     and plus-minus (FPM)
                         G:p s  decimals for gain and
                                level for format switch (fgain)              G:6 100.

Qpr_Drift    - logical - include drift parameters in the protocol.    [.false.]
                         They are signed so that the signal reduces
                         the gravimeter record. Beware the calibration
                         factor. The are formatted ready for tsfedit.
                         Fishing for the commands can be done with sed:
                         sed -n '/::/s/.*:: //p' prl-file
Qalways_To   - logical - If set .false., the tsfedit commands for the drift
                         terms will avoid those To-dates that end with
                         the time series and only write the From- and
                         eventually the Tref- dates.                   [.true.] 

                      RESPONSE COEFFICIENTS "TUNING"
_____________________  NEARLY DIURNAL CORE RESONANCE _________________________
Q_NDR       - logical - Effective tide potential includes NDR              [T]
                        If true, x_NDR and ZFacL2 apply:
Fre_NDR     - real*8  - Resonance beat period                         [435.d0]
Qua_NDR     - real*8  - Quality factor                                [2760d0]
ZFacL2      - cmplx*16- NDR-factor for Love L2                      [(6.0, 0)]
Core_model  - char*10 - NDR factorisation for
                        AR_STRAIN.  SP_STRAIN.  OT_EFFEC..        [SP_STRAIN ]

iun_prf     - integer - File unit from which a decimation filter 1)

                        will be input                                     [-1]

trg_prf     - char*64 - target label for the filter                      [' ']

dt_prf      - real*8  - The ratio of sampling interval prior to
                        decimation to the one after decimation         [0.0d0]
                        A value 0.0d0 means that the sampling rate
                        from the file is to be used.
                        Example: 60s decimated to 600s, dt_prf=0.1

o_prf       - char*8  - Retrieval option, 'B' and/or 'R'                 [' ']
                        'B' - Rewind before reading
                        'R' - Allow one rewind during search

thresh_prf  - real*8  - Threshold for printing, potential units.       
                        Negative numbers default to internal value
                        set in f_on_pot (0.5d-3 in potfcor.f)          [-1.d0]

qshow_prf   - logical - print the coefficients                       [.false.]

n_bessel_prf -integer - order and                                         [-1]
f_bessel_prf - real*8 - characteristic frequency of the analog
                        Bessel filter (GWR gravimeter). Off if n<0    [0.0174]     

1) A filter that has been used in the data preparation, decimating to a low-rate series like e.g.
  a 1h smp.interval. Read on here, esp. concerning cascaded filters and decimation schemes.

_____________________ AIR PRESSURE PERTURBATION ______________________________
Iun_SRT     - integer - file unit for *.TRS input (atmosph.job's tide
                        admittance table, output on file unit 12).        [92]
                        The atm.job's dump file is read from unit IUN_SRT-1
                        IUN_SRT < 0: No atm.corr.
Jaa_SRT     - integer - last argument number to check for association      [4]
Mxn_SRT     - integer - maximum degree of solid earth tide to charge       [2]
Zadmp(0:2)  - cmpx*16 - pressure admittance in record, tide band 0..2
Zadmt(0:3,2:3)  "     - a priori tide admittance, tide (order,degree)=(m,n)
PScal       - real*8  - units conversion [tide]/[airpress]           [0.01019]

_____________________ SPECTRAL ANALYSIS OF RESIDUAL __________________________
            Max.entropy part c.f. MXLPEF and Ls_Mem
            Method "Bartlett" or "Periodogram" is selected on the basis of
            data length versus spectral length:
            if ( Ls_Psp .lt. L_Sp/2 ) then Bartlett else Periodogram
            where L_Sp:=MIN0(actual_data_length, Lt_Pdg).
            Ls_Psp can be changed interactively.
            The feature is relevant only when urtapt is run with graphics on.

MxMiss_Psp  - integer - allow indicated number of missing samples          [0]
Wt_Psp      - char*4  - Window type (cf ~/sas/p/window.f)             ['KAIS']
Ls_Psp      - integer - Autocovariance length                            [256]
Wdp_Psp     - real*8  - Design parameter for Kaiser-window in resid.psp  [1.8]
Lpef_cfd    - integer - Length of PredErr filter for PSP confidence       [-1]

                        If <= 0 the sliding-average method is used.

Lt_Pdg      - integer - Length of input for periodogram. Data will be
                        windowed and zero-padded to achieve
                        length 2*L_Sp                                  [99999]
            - logical - Use a window for downweight. This window

                        is declared in *.INS file Downweight block         [F]

_____________________ DISPLAY _________________________________________________
q_shut_graphics - log - This option must be .true., the implementation of
                        refreshing a graphic window screen is faulty.

                        "TS" = time-series
DISPLAY     - char*12 - O - observed TS, also those under
                            "+" and "#"  tide band commands,
                            also multi-component file input,
                            and the columns of the normal matrix
                            ("observed" wave bands)
                        I - input TS under "+" tide band command.
                        X - input TS under "#" tide band command, after accomodation
                        # - input TS under "#" tide band command, after tsf-edit
                        U - updated TS after "+" tide band command.
                        E - updated TS after TSF-Edit

                        M - multi-component
                        B - wave bands

                        P - theor.potential (not LSQ fit)

                        T - theor.tide, LSQ fit, predicted and residual
                        t - theor.tide, LSQ fit, predicted

                        N - nontidal, LSQ-fit, predicted and residual
                        n - nontidal, LSQ-fit, predicted

                        R - all residuals (except filtered)
                        r - tidal    residual
                        L - least-squ.  "
                        F - filtered    "

                        W - weights
                        V - eigenvectors in model space
                        C - covariance matrix after eigenvalue editing
                        c - correlation matrix after eigenvalue editing
                        @ - Autocovariance of filtered and weighted residual
                            (before PSP)
                        A - all (but not W or #)                          ['A']

Display_gtg - char*8  - Controls the appearance of the GTG matrix as called
                        from subroutine BPROD (urtapu6.f).
                        +L - logarithmic display and automatic scaling
                        -L - linear display
                        +U - unity scaling under linear display
                        -U - raw scaling                                   [' ']
                        The subroutine's internal defaults are '+L-U'

            - logical - Time series display: Dots instead of line

npal1,npal2 - integer - From URTAP.PAL palette file, use line number
                        npal1 - for text screen,
                        npal2 - for graphic screen

List_TGP    - integer - list the strongest partial waves of potential

_____________________ TEST ____________________________________________________
Q_Frq_Tst   - logical - A pure cosinusiod replaces, injected after
                        whitening filter has been determined on real data   [F]
Frq_Tst     - real*8  - The frequency of that cosinusoid [cyc/d]         [1.d0]

_____________________ MODEL COVARIANCE ________________________________________
Q_transm    - logical - Process a priori model space covariance             [F]
                         Default = unity COVmm = implicit in SVD
IUN_transm  - integer - File unit for input of COVmm.                       [4]
                        C.f. Subr. Proc_Transm  in  URTAPU2.f

Thoughts about leap seconds

If input data is taken on the same UTC second (it's the time stamp that's relevant here), regardsless of intervening leap seconds,
then the data i not exactly regular in TT (Terrestrial Time). This is certainly the case in externally acquired data.
In order to generate tides exactly in synch, use the Namelist option  q_inject_leaps=.true.

How is that if external data is stamped with the Unix clock? Well, after a leap-second injection, the Unix clock will jump back.
So also here you should use     q_inject_leaps=.true.

There is only one circumstance when you should set .false.: If your data is surely equal-spaced. Therefore, the default setting is .true.

Example for a TSF EDIT block and instructions for processing and data weighting.
It's from TD/a/Allcamps/urtap-201502b.ins
(with a mistake, sorry). File
21 < d/

is multi-component (four columns) containing absolute gravity obs. and sigma,
and signal to adjust in regression (SCGrav obs. and drift)

 q_tsf_edit=.true., tsf_edit_name='OS '               <-
 tsfile_names='o/scg-cal-201502b '
 outlier_criterion=-5.0, iun_outliers=31
 q_allow_dc=.true., q_weights_dc=.true.
Tide Bands
#21,'BIN',-99999.0,'L:S|V',3,0,-1 = SCG
+21,'BIN',-99999.0,'L:S|D',3,0,-1, -1.0              
<-  that seems more appropriate for drift correction

End Tide Bands

Downweight Residual
!DT F=/0: 21,'BIN',-99999.,'L:A|S',3,0,0,0.
0 Normal

Downweight *
Q !DT F=/: 21,'BIN',-99999.,'L:A|S',3,0,0,0.
0 Normal

@21 +N L=
S|D A=-1.0                <- subtracts drift (what!! it's in Volt but applied to AG. Wrong!!!) 
CONT 31 O o/scg-cal-201502b.tse      

Tidal gravity:  Delta system
The response parameters are computed as follows:
The basic amplitude coefficients are tide potential elevation (TPE) [m] = m2/s2 / GE.
For obtaining response coefficients, the observed tide amplitudes are divided by the TPE.
In order to show delta-factors, we divide by RE/(N GE) . That brings us to an equatorial radius.
In the Dehant system the station's geocentric radius is used.
Although you might have chosen  IALTG=0 , you may ask for  DELTA_SYSTEM='DEHANT'  to get the traditional kind.
An illegal combination would be  DELTA_SYSTEM='DEHANT'  and  IALTG =/= 0


First, we give the most often used method for downweighting data using uncertainties from an input file.
The instruction file contains a Down-weight Block
The system searches at three stages for downweight information. Search names (Targets) are


There is an additional option to call the weighting procedure from a TSF EDIT block with a user-defined search name.
No, this option was dropped. It is more useful to call TSF_EDIT from the weight processor. We often
read an mc-file, where we get the weights from a data column. If we do things to the values we should
do equivalent things to the weights.

To enact down-weighting in the Normal equations, the following namelist parameters are significant (default values in [ ]):
Q_Downw = .true.   [.false.]

To enact down-weighting on the residual  
use_window_residual {'n','N'}    ['n']
use_downw_residual {'n','N','?'} ['?']  

and the block Downweight Residual must explicitly appear in the instruction file.

Description of Down-weight Block, file processing.
 Codewords and compulsory text are bold and underlined.

 Variable parameters appear in italic. Optional parameters between [ ].
 For full information, cf.

 Line#  Instruction              Explanation

   1    Downweight [target] [; comment]
                              Code from column 1, first eight chars spell exactly this way:
                              (DOWNWEIG or downweig or Downweig)
                              target - Instruction block identifier to match a search
keyword, any of NormalEq or Burg or Residual,
                                       or use * in order to match any search key.

         When Proc_Downw is called for acting on the NormalEq, treatment of DC-levels must distinguish
         between the data series Y(m) and the design matrix B(m,n)

   2  [ ?DT | !DT ][ {RDC|RWM} ][ DCA{+|-}] F[=operator]: lu,trg,recmrs,fmt,khms,itz,lddf,v [options]

       WFFI = weights-from-file instruction.
       Sampling interval, options, file unit, read and
filter options. 

                 File is processed with Proc_f_on_rec
                 (cf. ~hgs/sas/p/urtapu13.f),
                 weighting is carried out with Mult_ts_ts
                 (cf. ~hgs/sas/p/sasu061.f).

                         ?DT - check compatible sampling rate of data vs. weights
                         !DT - impose data's sampling interval on weight series.
                         RDC - remove DC-level in data (or residual) series before down-weighting.
                               Is often necessary !
                         RWM - remove the weighted mean from the data series
                        DCA+ - keep the DC-levels in the design matrix
                        DCA- - remove (except in a column that has been marked for carrying the
                               DC-level+  This is executed with a CALL Proc_Downw_DCC (n) )
                          F= - codeword
                    operator - string is passed to mult_ts_ts:
                           / - divide (else: multiply).
                           Q - take square root of input values.
                           2 - take square of input value.
                           0 - return weighted data with zero DC level.
                               Q and 2 are mutually exclusive, but 2 will win.

                 Additional instruction records might be needed due to filter options (due to lddf),
                 e.g. lddf = 0 => an additional line with a scaling multiplier is expected
                 immediately after the WFFI.

                 Additional options Q  - quiet
                                      DBG - debug

    2.1  [ TSFEDIT [ { O=un o filename] | F=filename | U=un } ] T=target [V=w] [+D] [+L] ]
                             - Optionally call tsf_edit, instructions from file.
T=target - mandatory.
               un , filename - Use an open file by name or unit number. No default for un.
                               Specify at least U=un ; the ins-file (usually 5) cannot be used!
            O=un o filename] - To open a file at this point, use the O= form.
                               If TSFEDIT is presented directly after Downweight (no file input), the
                         V=w   option can be used to create a constant array of value w (default: 1.0),
                               else w will be ignored. The edit commands will alter this array.
                          +D - Optionally, the ordinates deleted in the weights array are also
                               deleted in the data array. Useful only in the Residual case!
                          +L - By default the weight ordinates are deleted where the data ordinates
                               are missing. This can be inhibited with +L.

    2.2  [ SELECT { KEEP | REMOVE } DC :i[-j][,i[-j],...]] ]
                             - select columns in the design matrix for removing or
                               keeping their DC-level after weighting.
                               The column numbers are coded as ranges i-j
in the
                               list following the colon. The list must be coded
                               without blanks.
                               In the namelist, REMDC_Postweight should be set to '...'
                               (the default is 'AUTO').

    2.3  [ SPEC[ ]string ]   - string is sent to subroutine modify_weights. Withdrawn.

    3    0 Normal            - This section can be used for windowed or individual weighting.
                               It's maybe obsolete.
                               Instructions can be found in ~/sas/p/downws.f subroutine proc_downw .
                  `0 Normal´ - is the trivial case.
                         `0´ - is the number of records to follow that specify individual sample weights.
                    `Normal´ - is a codeword for using a window only in the building of
                               the normal matrix. Default is full-width rectangular, thus no effect.
                               Examples are given below.)

Long time ago we allowed window declarations. This option has been withdrawn (is discouraged). There are methods available anyway.
One is to synthesize a weight series or modify a copy of an existing file of measurement uncertainties.
Another is to create a file containing a constant value, e.g. unity, and modify it with a tsf-edit tse-file.
Still another method is to let subroutine proc_downw create a constant-valued weight array and let the user alter it
using a TSFEDIT command with no preceding file reading F=  (or F:) .


Here, the sigma-values from a GPS baseline LTU-file are used.
The two Downweight instructions below belong together:

Downweight Residual
F=/0: 21,'TOP] Rwd Y:92',-99999.,'(f12.6,t65,f10.0)',0,0,-1,0.
0 Normal

Downweight *
F=/: 21,'TOP] Rwd Y:92',-99999.,'(f12.6,t65,f10.0)',0,0,-1,0.
0 Normal

The first instruction is used in the RMS-residual stage, the second
in the other cases, which are `NormalEq´ and `Burg´.
They differ in the operator, `/0´ and `/´, respectively. The former asks
for a zero mean in the resulting weighted series. The `/´ asks for
division (inverse sigma = weight). The resulting residual RMS value should
be close to unity. It's the square root of the normalized Chi2.
In batch mode, the residual is exposed to weighting only if either
index('nN',use_window_residual).le.0  or
(namelist; default is  'N', and  use_downw_residual  has priority).

Downweight Residual
!DT F=/0: 26,'BIN',-99999.,'U',0,0,0,0.
0 Normal

because of the zero (shown in red)

Downweight *

F=/: 38,'BIN',-99999.0,'U',0,0,-1,0.


0 Normal

to modify file data from unit 38 synchronous with the observation array and use as weights.
The user has full control (responsibility) to convert the values into multiplicative, positive valued factors:
POKE 1.0
DEL From 2009 07 07 0 0 0 000 To 2015 08 30 23 40 00 000
Remember: these would be dividing weights. The example shortens the analysed time series cutting out the first part.
Since all filter operations will be done before the weights are applied and the columns contracted, any part of the series
could be blanked out.


            subroutine  proc_downw (Target,iuins,Display,

Refer to /sas/p/downws.f  Abridged information follows here:

Processes a "down-weight block" on unit iuins

Down-weights data, incl. application of a data window on the whole length
of the records.

Down-weighting of sections of bad data may be combined with a smooth
transition of weights (cosine-type).

Target     - char**  - to locate the instruction section. If empty, a lonely
                       "Downweig" is searched for.
                       Keep string short, six to eight chars.
IUINS      - integer -
Display    - char**  - If 'W' or 'A' and file data, call display_ts weights
B(M,N)     - real*8  - signal model
Y(M)       - real*8  - data
N          - integer - 0: no down-weight on B, down-weight on Y;
                          no window on B nor Y; overrules Block instruction.
                      <0: down-weight on Y, optional window on Y.
                      >0: down-weight (and window) on B and Y.
MU         - integer - signal length
EFFW       - real*8  - returned: factor to be applied on RMS(Y_returned)
                       to achieve unity gain of weighting process.
T0,DT      - real*8  -
WRK(NWRK)  - real*8  - work space, min. length: slightly larger than mu



  Proc_downw_tsf_edit (iun)     - specifies a file unit, a file with a
                                  Tsf_Edit block for
                                  processing of the weight series

  Proc_downw_rwmean (qq)        - .true.: remove DC in Y(.) by
                                  Weighted_Mean (sas/p/wmeans.f)

  Proc_downw_rwmy (rwmc,drwmc)  - returns the weighted mean of Y(.)
                                  and the sum of weights WRK(.)

  Proc_downw_dcc (n)            - defines the DC column in B(.,.)

Description of Down-Weight Block: refer to ~/sas/p/downws.f.  


Unless careful consideration, do not request windowing with HANN, KAIS, ...
Their mainlobes are simply too wide. Unless you accept to sacrifice tide
wavegroup resolution. If you have a long record, try first to separate
K1 S1 P1 with the rectangular window. Then you might try to combine
K1+S1 in one band and use the HANN window.
   The quarter-taper window HANQ might provide an acceptable trade-off.
It maintains the rectangular window spectrum at small frequency offset,
whereas sidelobes approach those of the HANN window at large frequency offset.
   Do not consider tapered windows when you have to down-weight outliers.
These weights upset the resolution capability much more than what could be
counterbalanced with a window. Only extreme outliers should be down-weighted.
   More explanations are found in /sas/p/downws.f

Examples for Downweight using explicit weights:

Downweight * ; in all instances the following weights are applied:
6 T=9                            ; cosine taper = 9
89 07 10 14 00 00  227 0.3
89 09 23 20 00 00   81 0.2
89 10 20 14 00 00   73 0.4
89 11 05 17 00 00   32 0.05
89 11 24 02 00 00   41 0.2
89 12 02 20 00 00   50 0.2

Downweight * ; does not do weights in the residual RMS section
7  Window=HANQ in Normal Equations        ; HANQ-window the whole lot of data
89 07 09 08 00 00  234 0.25
89 08 22 11 00 00   41 0.3
89 09 24 17 00 00   53 0.1
89 10 10 05 00 00   35 0.2
89 10 21 17 00 00   55 0.2
89 11 06 02 00 00   32 0.05
89 11 30 11 00 00   41 0.3

      subroutine proc_f_on_rec (crec,iuins,iun,w,nyo,n,t,dt,v)
/sas/p/urtapu13.f   might be newer. Unmaintained commentary follows:

Read a time series from a file (unit IUN). process with a filter.
Take parameters for reading from string CREC.
Read filter control instruction string from file unit IUINS.
Process time series with the filter.

CREC    - char*80 - Command line; replaced by filter line on return
IUN     - integer - log.unit, returned
W(nyo)  - real*8  - Data array, returned, original dimension NYO
N       - integer - Length, returned
T,DT    - real*8
V       - real*8  - c.f. below, returned to calling program

Command line:
IUN     - integer - file unit
TRG     - char*32 - target string for file positioning in the case of ASCII
                    files or file type codes e.g. 'BIN' | 'BRX'.
REC_MRS - real*8  - missing record symbol
FMT     - char*64 - format string | 'U' case 'BIN'ary files
KHMS    - integer - time record check (1..3)
ITZ     - integer - time zone (+1 for CET)
LDDF    - integer - filter option
V       - real*8  - constant. Needed for '+' (LINCOMB) processing
                    Default = -1.

        FMT or TRG may contain additional instructions like
               'Rwd'  - rewind file before reading
               'Y:92' - set year 1992 for decimal year date format.
               This is passed to Read_FuF / Read_Fm_D.

               or -2: Process with scaling.
                        The command line following
                        contains the scale factor.
                        Optionally, the scale factor may be used for
                        converting air pressure into input-data units;
                        this is effectuated if the line contains 'PCAL'.
                        Recommended if air pressure results will be
                        combined into the present solution.
              -2 or -3: Subtract DC-level
            -1 or < -3: Take series as it is. No command line follows.
                   > 0: Next command line describes filter:
                        Scale factor, NFM,NFP,SYMM.
                        Follow: NFP-NFM+1 filter coefficients.


           subroutine  proc_rslist_add (nie)

Reads *.INS file and requests printout of minor tides from RSLIST

Control block code word = 'Show Also...' starting in the first column
Data block:
(1)   N = number of lines following, read (4,*)
(2..) Q0...Q7 SYMBOL               format (8i2,1x,a4)

where Q0...Q7 are the astronomic argument numbers and SYMBOL the Darwin
symbol of the requested tide.

N.B.: As usual, Q0 = spherical harmonic degree of the tide,
      Q1...Q6 the argument numbers for tau,s,h,p,N',ps,
      Q7 the cosine/sine/sign-selector.

            subroutine  proc_del_Xtalk (iun,a,n)


Processes named Instruction block "Delete X-talk" on file unit 4.

Deletes off-diagonal elements in normal matrix. (Doubtful value.)

Instruction block:

Delete X-talk                - codeword, spelt exactly this way
N                            - number of records following
i1 i2 j1 j2                  - do-loop limits for off-diagonal elements
                               N such records

                             Data is read using READ (4,*) ...

(Comment: This routine is nearly deseased. It gives strange results and
the name might raise hopes it can't fulfill. We should delete the option.)

            subroutine  proc_transm (iun,v,n)

Processes *.ins file records and sets model covariance matrix

v(n,n) work array,
iun    file unit for value substitution:

File should contain the following:
Column 1----\
Line 1      Model Covariance              = identifier string
Line 2      n                             n=number of records (2 lines):
Record 1.1  i1 i2 flag v1 v2 k
Record 1.2  d(i1) d(i1+1) ... d(i2)  if  flag=F, real SOLP(i1..i2)
alt.   1.2  d(i1) d(i1+2) ... d(i2)  if  flag=T, complex SOLP(i1),SOLP(i1+1)
                                                     ... SOLP(i2),SOLP(i2+1)
       1.3  i j                      array index (i<j) for sign reversal
...  1.2+k  i j

Then, diagonal element = d(i)
off-diagonal           = root(d(i)*d(j))*v2*(v1**iabs(i-j))
for  i1 <= i <= i2  and  i1 <= j <= i2.


(m = model parameters (space), d = observed data, B = signal model array)
(Text = user entries at prompts )

Phase 1   Input data.
          Read d, d1,d2,... and compose  d <- d + d1 + d2 + ...;
          If Multi-channel input data: prompt to select channel (by name).
          Assemble B, possibly adding a detrending signal and a constant.
          Dump model array and observations.

Phase 2   Get dump;
          Filter columns of B and d, using  DFFF <- DFF*F;
          Optionally apply weights and/or taper (Process "Downweight" block)

          Prompt user for solution options
          +S - store previous solution parameters: mold
          R  - reset mold
          +U - solve Delta_m = mnew - mold  from Delta_d (the previous
          +T - transform w.r.t  m mT-covariance

          make Normal equations;
          optionally transform w.r.t. m-covariances (using Jocobi-routine).

          Solve inverse system by Jacobi-eigenvalue routine.
          Prompt user do remove small eigenvalues (singular values).
          Show model resolution and VarCovar in transformed space {m'} or
          original space {m}.

          Solve parameters (and transform  m' -> if necessary);
          Print results table (actually postponed to after Power Spectrum;
          we need the noise-RMS from the residual analysis).

Phase 3   LSQ residual. Optionally display and...
          optionally output  d - B mtides       "Tidal residual"
               "       "     d - B mall         "Residual"
               "       "     DFFF*(d - B mall)  "Filtered residual"
          Burg prediction error filter determination (PEF), using residual.

          Spectral analysis of filtered residual:
          Prompt user to modify spectrum length and window parameters.
          (Enter any character, ICHAR < 64, to stay prompted after spectrum
          display with initial parameter set.)
          Bartlett Power Spectrum or Periodogram if data too short; indicated
          on screen plot only by y-axis label: Power or Amplitude, resp'ly.

          Maximum entropy spectrum:
          Show power spectrum of  PEF*DFF; prompt user to select appropriate
          order of PEF.

Phase 4   GO | MC-next | STOP
          Prompt user to decide: GO = rerun from phase 2, or MC-next or STOP.

MC-next   Multi-channel input data: Select next channel and goto Phase 2.1

GO        PEF*DFF -> DFFF, goto Phase 2.

STOP      Write *.evs-file, *.trs-file and stop.
          Program can be reissued from Phase 1 or 2 (Qget_dump = .false. or
          .true., respectively)


Q_transm=.false.  or m-COvariances small, deleted eigenvalues:
       confidence limit increases, solved parameter obtains uncertain value.
       Remedy: introduce m-covariances and solve with Q_transm=.true.
       or:     withdraw the parameter
       or:     subtract the effect from the observations if an admittance value
               can be inferred with decent accuracy from elsewhere.
       or:     don't delete the eigenvalue

Solve with Q_transm only if reliable model covariances can be supplied.
A good estimate of m should have been achieved before solution option
+U is issued. Weak parameters might obtain large confidence limits; however,
pushing the system to solve a weak parameter "in the mid" of the confidence
limit is supposed to results in low bias of the other parameters.
  For example, an L2 admittance should be near S2 and M2. Using -U,

m(M2) = 0.31 +- 0.01
m(S2) = 0.32 +- 0.02 and
m(L2) = 0.01 +- 0.5  is not wrong. Interpreting this on Delta-m, assuming
                     m(L2) ~= m(S2) would be more satisfying though:

m(M2) = 0.31 +- 0.01
m(S2) = 0.31 +- 0.02
m(L2) = 0.31 +- 0.5  using +U. The first solution must be "good" for this

Example for TSF EDIT:

           * Open file block
        24 < ~/gps/ap/onsa93.met.ts
        23 < ~/gps/ap/wett93.met.ts

        Tide Bands
        #23,'BIN]',-9999.9d0,' ',0,0,-1 E:APRD] = aprO
        End Tide Bands

        24,'BIN]',-9999.9,' ',0,0,'L',-1.d0
        0.0,0.5,0.0                  Wettzell pressure 2 times as efficient
        DUMP 29
        SCALE 0.001  to get responses in mm/hPa

Tide bands section calls reading procedure to input file #23 and TSF Edit
procedure to operate on, tag APRD.
TSF EDIT section tag APRD requests reading of a second file, subtract the
data with a scale factor of 0.5, remove DC-level, and dump the result on
file unit 29.

URTAP can perform Empirical Orthogonal Function analysis. It can also
treat mixed problems of fitting deterministic and empirical functions

An example is the fitting of several residual series of neighbour GPS
sites to one set of observations. By this we hope to delete common mode
signal patterns that persist in a region or are imposed by a subset of
the orbit/regional network constraining sites. The common mode is represented
by the largest eigenvalue and the associated eigenvector.

In the wave group section we load data

##iun  open file statement
# proc_f_on_rec statement -> EOF
many more

We delete the 2'nd, 3'rd etc. eigenvalues that belong to this segment of
the equations:

                 eof_grades(2)=1.d0, eof_grades(3)=2.d0, ...
or, if you dare,
                 use_eigenvalues='D n G '

where n is the number of EOF's for which the eigenvalues are to be deleted,
usually the number of input series minus one.

In the latter case we suppose that we know that eigenvalue #2 belongs to
the EOF-part and is the second largest in magnitude. In the former case we
can have a look at the associated eigenvectors and delete them by plausi-
bility. There is a little structure-function analysis that may assist
in finding nearly-white-noise time series. That part of the program is
certainly worth further development.

The Display option V will show the model-space eigenvectors. They are
useful to inpsect the degree of participation in the common mode and
the modes that are orthogonal to the common mode.

The Display option c will show an array of correlations between the
data-space eigenvectors after EV-deletion. The diagonal is identical
to the sequence of resolution values.

The auto_select_eu method may fail to identify sufficiently many eigenvectors.
The criterion is to find vectors with a projection onto the D-space
exceeding 0.5 (The D-space is the space created from the nexternal time-series.)
If we instead use 0.5-m, the number of vectors will increase. The parameter
m can be specified through the namlist: auto_detect_eu_margin.

Decimation filter / Pre-filter


You have decimated the data series with a cascaded low-pass filter + decimation, e.g.
from 1s to 600s in two steps. In deci.tse for instance:

FILTER D:WD KAIS 40 2.1 0 0.05 1.0 0.0
FILTER D:WD KAIS 58 2.1 0 0.0333 1.0 0.0

For estimating unbiased tide coefficients you must input this filter using the
namelist options *_prf in the form of its coefficient series.
Generate the filter BIN-file using

cd ~/TD; tslist _3000 -B -rs1. -Eplot/prf.tse,600 -N | m

with  prf.tse re-engineered from deci.tse (changes appear in fat)

POKE 1.0 At #0
FILTER +C D:WD KAIS 40 2.1 0 0.05 1.0 0.0
FILTER +C D:WD KAIS 58 2.1 0 0.0333 1.0 0.0
OUTFS X> +S U=31 < prf.fs] From #0 To #58

The "#58" comes from the length of the second filter.
Insert a DEBUG command before OUTFS if you are uncertain.
The second DECIMATE command must be deleted.

Concerning tslist, the 3000 samples are sufficient to accommodate the filter. Check that you do get a central section of zeros.
The -r option must specify the sampling rate of the input data.
The -B option is only a nuisance.

We have an alternative to use a long low-pass filter in the spectrum domain for decimation from 1s to 1h in one quick step:
You must provide one day's 1s-file with two hours pre- and appended. That's easily done with ncd2ts,
learn from section (3) in
    ncd2ts -h
With the three ncd-files in shell variable $ncds and the centre date in $ymd 
    ncd2ts +% -E deci.tse,SPF1H -022. -UEn100800 -o tmp/G1$ymd-1h.ts % -c Grav-1 $ncds

Back from where you came

Drift model output
, tsfedit commands


urtapt @ urtap-bigtidesww.ins

sed -n '/::/s/.*:: //p' oa/g100201-130515-1h-lapww.prl
ADD  -3.25393464D+01 From 2011 02 24 01 00 00 000  To 2013 05 16 19 59 59 000
SLO  -1.14597635D-02 [/h] Tref 2012 04 05 10 30 00 000  From 2011 02 24 01 00 00 000
SLO   1.38059881D-02 [/h] Tref 2011 09 24 21 30 00 000  From 2010 02 02 00 00 00 000

Place these lines into a  tse-file, wrap with TSF EDIT TARGET and END

tslist G1_xyz-10s.ts -E tse-file,TARGET ...

Back from where you came