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).
o
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/tide-phasors
SCG/plot-vcv
(binary vcv-matrix)
SCG/plot-vcvf
(ASCII)
SCG/plot-solrms
(2) There is a help file HOW.TO
GPS in which urtap applications appear.
Multiple time-series fitting may involve an
Empirical Orthogonal Function
approach.
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.
o
The program produces a couple of output files, time series
containing
(1) the tidal prediction and the tidal residualSome of these can be deselected using namelist options. In addition, aowhole suite of
(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.
Instruction file structure:
===================
1 Open file block
2 Namelist ¶m
3 INSTRUCTION BLOCKS
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
o
TSF EDIT what
codeword and tag
END
o
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
urtape2.f
Character strings that are examined for environment
variables:
================================================
- All file names in the open-file block
- Namelist: fmt, tsf_edit_name,
tsf_late_edit, tsf_edit_wr_name, tsfile_names,
core_model,
name_site, s_site, s_idate
- Instructions for wavegroups and signals
- However, not (yet): Instructions for weighting
Syntax:
${VAR}
${VAR:[default]}
_________________________________________________________________________________
o
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)
====================================================================
Syntax:
Every line that identifies a signal (although not always shown) must end in symbol assignment (one to four characters)
- Underlined characters are fixed vocabulary; they oannot be modified.
- Italics are to be replaced with values or text strings.
- [ ] indicates optional code.
- ... expected content that will be explained further.
- All declarations must start in column 1. Unidentifyable code will be ignored.
- Use a `;´ in column 1 to comment your coding.
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.
(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
(...)
| <...)
TIDE WAVE GROUPS or INCLUDED TIDE WAVE SUBGROUPS
(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
The logical unit for reading can be changed at any
place within the block:
i #iu
[o filename]
If only the unit number is given, the file must have
been
opened in the open-file block at the top of the ins file.
When the EOF on iu is
reached, reading continues
from the original instruction
file.
OBS! Cannot be nested further.
____________________________________________________________________________
Select signals
for a customized (hand-selected) predictions:
Bracket signal declarations with BEGIN_HS ... END
commands:
BEGIN_HS [n1 [n2 ...]]
selectors
BEGIN_HS [n3 [n4 ...]]
...
ENDBEGIN_HS [set-number [set-number...]]
selectors
BEGIN_HS [set-number [set-number...]]
...
END
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 (e.g. tsfedit-generated) or observed time
series, multiplied by the estimated admittance).
Presently this does not include commands
P and § (parabolas); these commands must
appear last before END
(perhaps not even Slopes???)
In an upgrade of urtapt.f , the maximum number of
hand-selected sets is scalable
( include file urtapuc.fp , parameter msxhs ).
Concerning Hand-selected sets in CONDUMPS:
Since we wouldn't like to output unused memory in the various dump
operations, an option is available:
q_hs_as_needed
(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
BEGIN_HS ... END.
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 (!), inverse selection is also possible:
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.
A side effect of detrending in form of a DC-signal may suggest to
set q_allow_dc = .true.;
it might even be decided internally owing to complex circumstances
(search for the string "DECIDE_DC" in the source
code). However, as another side effect, the design matrix
may become singular since a DC-column encompassing
the whole stretch may already have been introduced. In order to
prevent this situation, the automatic addition
can be prevented by setting q_allow_autodcband = .false. (is
.true. by default).
Since detrending may be selected as automatic,
you can decide where you include the trend
using the opt_hand_select_trend (namelist)
option:
opt_hand_select_trend = 'string'A partial residual can also be produced. using the opt_hand_select_resid (namelist) option:
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'
opt_hand_select_residThe hand-selected output time series are named
If the nk'th position of the option string contains a 'T', a partial residual time series will be
produced (observation minus prediction).
( rootname is assigned by
tsfile_names='rootname' in the namelist.)
____________________________________________________*______________________
\ o *
Coding TIDE WAVE
GROUPS:
______*______
Identifier = (
(argument selector) = tide-symbol [W=#iun]
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.
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
W=#iun - Write this band
to file unit iun. BIN mc-file with label tide-symbol
to
which _REAL resp. _IMAG will be appended. File
must have been opened.
Suitable for diagnostics, the other tide band commands
won't need this capability.
Either
files are present or output commands can be added when tsf-edit is
invoked.
Coding INCLUDED TIDE WAVE SUBGROUPS:
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
____________________________________________________________________________
Frequency in cyc/dIdentifier = (~
(~ frequency) = symbol
(~ 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 sheet.
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.
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
Equivalence cases for TARGET's extension (last two characters):
uniqueness the first extra letter should be `_´
so that no ending of a four-letter-only symbol
would match a two-letter target extension.
_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
Example:
= #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 and read
at "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
${SITE}_^
____________________________________________________________________________
Using TSF EDIT as a SIGNAL GENERATOR
Identifier = t
t #iutse target = symbol [+LET]
(You can code iutse as a string #nn
where nn is a two-digit number.
You can specify it without the `#´ too).
The file containing the signal generating commands either
must have been opened in the file-open block at the top of the
instruction file or, in particular if
you're going to exhaust the number of files that can be open
simultaneously, within the Tide Bands block:
##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.
Option +LET allows the case of a target that
cannot be located in the file. If +LET is not
specified,
the program will stop instead.
________________________________________________________________________________
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,lddf [-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
symbol.
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]
f(nfm),f(nfm+1),...,f(nfp)
where:
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=namemust 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
name - 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
Example:
#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:
##openf-string
#iuninx,...
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:
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, result sheet). An
"explain"-feature is under 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/jpl_geoc.mc
# 23 L:X_VAL D:0,0,0 R:.] = GEOX
0.001
##23 # /geo/hgs/geoc/jpl/jpl_geoc.mc
# 23 L:Y_VAL D:0,0,0 R:.] = GEOY
0.001
##23 # /geo/hgs/geoc/jpl/jpl_geoc.mc
# 23 L:Z_VAL D:0,0,0 R:.] = GEOZ
0.001
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.
The same time series can be read in more than one time, each
associated with a different symbol:
open file block:
25 ^ d/BPRIOS090615.ra-1h.ts
Tide bands
...
#25,'BIN',-99999.0,'U',0,0,-1 = RIBP
#25,'BIN',-99999.0,'U',0,0,-1 = OSBP
Then, at the stage between filtering and gap
contraction the time series can be made disjunct using namelist
iun_tsf_late_edit=32
col_edit(1)='RIBP',
col_edit_trg(1)='RIBP]'
col_edit(2)='OSBP', col_edit_trg(2)='OSBP]'
with the instruction file
32 B d/RIOS.tse
containing
TSF EDIT RIBP
REPDC From 2015 06 25 00 00 000 ->S
POKE % From #0 To 2015 06 24 23 45 000
END
TSF EDIT OSBP
REPDC From #0 To 2015 05 24 23 45 000 ->S
POKE % From 2015 06 25 00 00 000
END
Thus in the solution, each part - in this case the bottom
pressure series from tide gauges Ringhals and Onsala -
will obtain each their admittance coefficients
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
TSF EDIT LEKS
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
END
For early editing, the Namelist ¶m would contain
tsf_edit_name='LEKS]', q_tsf_edit_y=.true.
Additionally a late editing option has been made
available.
For backward compatibility the parameter q_tsf_edit
is recognized as the early-editing alternative.
For late editing,
tsf_late_edit='LEKS]',
q_tsf_edit_late=.true. [sic!]
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 log (unit 6) reads
<Main-->>>
No outliers detected.
File unit numbers:
In the instruction file at TSF EDIT CLEAN /
CONT 31
input unit number 31 is used merely because we know the
unit is available.
In iun_outliers=31 , unit number 31
implies that 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.
The search for outliers comprises the full range of the
(filtered) residual by default. However, in some cases
end-transients appear. To avoid scanning the tail of the series,
specify a (small) positive number in
namelist ¶m ,
e.g. n_noutliers_endcut=48
To write
outliers, the namelist
¶m would contain
tsfile_names='d/LEKS.ra '
iun_outliers=31, outlier_criterion=-5.0
Determining outliers in the filtered series (.wr.ts) implies that
the sample to be marked for delete must be located
some sampling intervals earlier; one sample because of the onset
of the filter, and perhaps one more due to a run-up
of the signal towards the extreme value.
This can be modified using e.g.
q_outliers_filtered=.true.,
idel_outliers_shift=-2
The default is -1
A strategy may consist of repeated outlier detection from the
unfiltered residual until no more are found;
copying the tse-file to a safe place; switching to the
filtered residual and find a criterion that yields parsimonious
deletion.
Filtered residual (.wr.ts):
In the preparation of the filtered residual there are two filter
stages, first the order-1 difference filter, then a PEF filter.
One-sample gaps can be closed by linear interpolation before the
first filter stage.
Before the second stage, a user devised repair can be applied.
Here e.g. linear interpolation of gaps max. 6 long.
q_close_gaps_wr=.true.,
opt_close_gaps_wr='L<=6'
Other options:
' <=#n'
- insert DC-level for gaps max. #n long
' ' - (empty)
insert DC-level for gaps of any length
Default is 'SKIP'
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 or 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.
________________________________________________________________________________
(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 Bandsdevises a partial residual which will be composed of
(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
B
= #22 (Site events file, look for) -> VANE_U
E
End Tide Bands
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
tsfile_name.ph.ts - 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
qouts_i_sel=.true.
________________________________________________________________________________
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] ] 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
Operation:
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
(2,0,*)
(2,1,*)
(2,2,*)
* Switch to SRT
(1:2,1,1,*) = S1
(2,2,2,*) = S2
#24,...
+22,...
+23,...
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
component.
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:
============
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 or ASCII
(qout_vcv_bin=.true.|.false.)
15 alternative
with
qout_vcv_bin=.false.
*.vcv-file ASCII
19 -
Solution array and design
matrix
*.gmx-file BIN
IUMCR - 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.
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} '
URTAP_LDATE - if namelist ¶m specifies
LDATE(1)=-1, $URTAP_LDATE will be used.
Specify year,month,day
URTAP_LDATA - if defined in environment, will
overcore namelist ¶m l_data
URTAP_NEXTRA - if namelist ¶m 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,
${name:[default]}
tsf_edit_nameE.g. tsfile_names = 'o/g${URTAPFD}-${URTAPDT:[1h]} '
tsf_edit_wr_name
tsf_late_edit
tsf_edit_weights
tsfile_names
s_idate
fmt
________________________________________________________________________________
File urtapuc.fp: Array size for signal bands
parameter (mxmcd=200,
maxhs=200)
MXMCD - Column dimension
(designates complex variables)
MAXHS - Maximum
number of Hand-selected sets
File tide_spectrum.fp: Dimensions of (Tamura's) harmonic tide
potential
parameter (maa=9,naa=1300)
MAA - Number
of astronomical arguments
NAA - Number
of partial waves
File wbands.fp: Signal band structure in
COMMON /C{C|F|W}BAND/
COMMON /CWBAND{C|P}/
parameter (ndft=1000)
NDFT - A
generous length for 13 arrays.
Main program:
parameter (nyd=1000000,
nbands=mxmcd, mxspef=500)
NYD
- max. time series length.
NBANDS - number of
signal components
(2
for
every
tidal
wave group).
MXSPEF - maximum order of PEFilters.
parameter
(ntotal=nyd*nbands+nyd)
NTOTAL - total size of design
matrix, columns * rows
parameter
(nbanda=nbands+10)
NBANDA - granting a safety margin (?)
in arrays {A|D|P|F|I|T}DET
parameter
(nvcv=nbands*nbands+nbands+nyd)
NVCV - Size of
variance-covariance matrix VCV with extra space
________________________________________________________________________________
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 applying
a
window; overrides use_window_residual
['?']
(Was introduced since use of windows is
discouraged.)
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, tailored for AG
multi-campaign analysis: 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
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, q_outliers_ddiff)
iun_condump - integer - if > 0 dump to this file unit
(must be opened in
open file-block
[-1]
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)
q_conmrs_weights
- logical - condense (and output) the weights. Automatically
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 ... )
iun_sbsel_out
- integer - File unit for dump of selected columns of the signal
array
right after filtering and before
condensation [-1]
k_sbsel_out(200)
- 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,...]
opt_sbsel_out
- 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 k_sbsel_out.
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; external
observations may be specified in
the `Tide Bands´
block.
['TGP']
Demodulation_Tide
- 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')
Eff_Average_Time
- 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_allow_autodcband
- logical - .false.: Prevent automatic addition of a DC-signal
in the design matrix. The normal matrix might
become singular otherwise.
[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]
Q_Tref_Slopes_Start
- 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
[F]
2. subroutine conmrs
[F]
Col_Edit(100) - char*4 - Up to 100 symbols
designating signal columns
(not: tides) for tsf_edit with command sections
in the file on unit
iun_tsf_late_edit.
Target
strings must be given in the respective
Col_Edit_Trg - char*32 - (dimension 100). Defaults are
'XXXX' and empty.
This editing occurs between filtering and
gap contraction. The subroutine used is
Edit_Signal_Bands (urtapu23.f)
_____________________ 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
[/home/hgs/tap/d/ttutcf.dat]
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]$']
qdebug_tidebands - log - Debug option for file reading in the
signal
model assembly
[F]
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']
Q_Del_Miss_Obs
- 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.
Opt_missalign
- 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]
Q_Burg_nice_slice
- 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]
REMDC_Postweight
- 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
[F]
2. subroutine conmrs
[F]
Q_Analyse_Slopes
- 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_________________
iun_outliers - integer - file unit to which
outlier records are written.
outlier_criterion
- 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]
q_outliers_filtered
- logical
- .true.: Scan the filtered and weighted residual
.false.: Scan the raw
residual
[T]
q_outliers_ddiff
- logical
- .true.: Scan the raw residual after double-
differencing
.false.: Scan the raw
residual
(q_outliers_filtered must be
.false.)
[F]
q_outliers_bypos
- logical -
.true.: write TSF EDIT commands with position
indexes
.false.: ... with date
records
[F]
n_outliers_cutend
- integer - No outlier detection in the tail of the
series [0]
qtalk_woutliers
- 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
where 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.
Q_TSF_Edit_Late
- 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
Q_TSF_Edit_wr -logical - At the stage when the filtered
and Weighted Residual
is inspected for the RMS
(significant in the
result sheet), Tsf_Edit can be
called to e.g.
delete samples that have been
deleted in a
down-weight block. See the example.
[F]
IUN_TSF_Edit_wr
- integer - The file unit number. Must have been
opened in the
Open-block. Default 0 refers to the same file
as the TSF_Edit_File
[0]
TSF_Edit_wr_Name
- char*64 - The target string for the edit
block
[' ']
_____________________ EIGENVALUE ANALYSIS ________________________________________
qfail_on_evcheck
- 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.
[F]
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
eigenvalues/-vectors.
Dimension eof_grades = nbands [50]
[nbands*0.0d0]
q_user_select_eu
- 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]
q_auto_select_eu
-
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.
auto_select_eu_margin
- 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]
Q_ra_with_outliers
- logical - The residual and the filtered residual will
contain the outliers (i.e. the observations
before they are processed by tsf_edit)
[.false.]
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
['..........']
opt_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
[.false.]
qout_vcv_bin - logical .false. : ASCII output else BIN
[.true.]
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
['ev-dspace.mc']
_____________________ 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.]
Fn_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] NEW IN NAMELIST !
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]
Slope_time_base
- char*8 - Representation of slopes in result sheet,
e.g. units per year.
Possible values: 'year'
'day' 'cty'
cty =
centuries
['year']
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
]
_____________________ PRE-FILTER
= ANTI-ALIASING IN TIME SERIES PREPARATION
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]
Q_window_residual
- 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'
Q_Display_Ts_Dot
- 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
obs.data, 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
¶m
q_tsf_edit=.true., tsf_edit_name='OS '
<-
tsfile_names='o/scg-cal-201502b '
outlier_criterion=-5.0, iun_outliers=31
q_downw=.true.
q_allow_dc=.true., q_weights_dc=.true.
...
&end
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.
1.0
0 Normal
Downweight *
Q !DT F=/: 21,'BIN',-99999.,'L:A|S',3,0,0,0.
1.0
0 Normal
TSF EDIT OS
@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
END
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
Down-weighting
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
NormalEq
Burg
Residual
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
.AND.
use_window_residual .not.in. {'n','N'} ['n']
.OR.
use_downw_residual .not.in. {'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. /sas/p/downws.f.
__________________________________________________________________________________
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).
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').
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:) .
Examples:
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
index('nN',use_downw_residual).le.0
(namelist; default is 'N', and use_downw_residual has priority).
Downweight Residual
!DT F=/0: 26,'BIN',-99999.,'U',0,0,0,0.
774.4
0 Normal
because of the zero (shown in red)
Downweight *
F=/: 38,'BIN',-99999.0,'U',0,0,-1,0.
TSFEDIT U=30 T=DW
0 Normal
to modify file data from unit 38 synchronous with the observation array and use as weights.TSF EDIT DW
The user has full control (responsibility) to convert the values into multiplicative, positive valued factors:
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,
b,y,m,n,mu,effw,t0,dt,wrk,nwrk)
=======================================================
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
____________________________________________________________________________
ENTRIES:
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.
Experience
==========
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.
LDDF
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.
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
Operation:
(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
residual)
+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' -> 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)
Experiences
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
purpose.
Example for TSF EDIT:
---------------------
* Open file block
...
24 <
~/gps/ap/onsa93.met.ts
23 <
~/gps/ap/wett93.met.ts
...
Q
Tide Bands
...
#23,'BIN]',-9999.9d0,' ',0,0,-1 E:APRD] = aprO
...
End Tide Bands
TSF EDIT APRD
24,'BIN]',-9999.9,'
',0,0,'L',-1.d0
0.0,0.5,0.0
Wettzell
pressure
2
times
as
efficient
REMDC
DUMP 29
SCALE 0.001
to get responses in mm/hPa
END
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.
EOF-analysis
------------
URTAP can perform Empirical Orthogonal Function analysis. It
can also
treat mixed problems of fitting deterministic and empirical
functions
simultaneously.
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:
q_user_select_eu=.true.
or
eof_grades(2)=1.d0, eof_grades(3)=2.d0, ...
or, if you dare,
q_auto_select_eu=.true.
anyroad,
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:
TSF EDIT 600
REPAIR L<=10
FILTER D:WD KAIS 40 2.1 0 0.05 1.0 0.0
DECIMATE 20 ${MOVE} 0
FILTER D:WD KAIS 58 2.1 0 0.0333 1.0 0.0
DECIMATE 30 ${MOVEKWD30} 27
END
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)
The "#58" comes from the length of the second filter.TSF EDIT 600
POKE 1.0 At #0
FILTER +C D:WD KAIS 40 2.1 0 0.05 1.0 0.0
DECIMATE 20 ${MOVE} 0
FILTER +C D:WD KAIS 58 2.1 0 0.0333 1.0 0.0
; DECIMATE 30 ${MOVEKWD30} 27
OUTFS X> +S U=31 < prf.fs] From #0 To #58
END
Example to restrict the analysis to a late part of the
observations Y (urtap-openend-esplc-cmhy.ins).
In the open-file block:
30 R urtap-cmhy-exp.tseIn the Namelist:
q_tsf_edit_wr=.true., iun_tsf_edit_wr=30, tsf_edit_wr_name='RW${CUTTO}'
Towards the end of the ins-file:
Downweight Residual
TSFEDIT U=30 T=DW${CUTTO} +D
0 Normal
Downweight *
TSFEDIT U=30 T=DW${CUTTO}
0 Normal
TSF EDIT DWgbalt
POKE 1.0
DEL From 2016 12 31 23 40 00 000
END
Content of urtap-cmhy-exp.tse to set the window
for the analysis:
TSF EDIT DWgnwsSetting the environment variable CUTTO:
POKE 1.0
DEL From 2009 07 07 00 00 00 000 To 2014 04 01 00 20 00 000
${DEL4GRWA:[; ]} From 2014 04 01 00 40 00 000 To 2015 08 31 14 20 00
DEL From 2018 12 31 23 40 00 000
END
TSF EDIT RWgnws
DEL From 2009 07 07 00 00 00 000 To 2014 04 01 00 20 00 000
${DEL4GRWA:[; ]} From 2014 04 01 00 40 00 000 To 2015 08 31 14 20 00
DEL From 2018 12 31 23 40 00 000
END
setenv CUTTO gnwsNow, the filtered residual will be cropped in the same way as the input data.
Drift model output, tsfedit commands
------------------------------------
qpr_drift=.true.
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 ...
.bye