Tide Analysis Package o
o
This program was originally created for analysis of ground water
records.
Resulting parameters are computed in terms of tide 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.
Quite generally, a number of tide problems can be solved
(CAUSE='TGP', CAUSE='TGP|SRT', CAUSE='TGP+SRT'). The program does
not compute the traditional response parameters (gravity: delta),
however.
What you get is an admittance between the observed quantity and
the tide
potential elevation [m].
o
Air pressure tides can be analyzed (CAUSE='SRT') operating URTAP
on a baro-
metric record. The results of the "Atmosphere job" can be processed
together
with the ground water analysis (CAUSE='TGP+SRT'), given odmittance
para-
meters for <ground water>/<air pressure> [m/hPa] and preliminary
admittances
for <ground water>/<tide potential> [m/m].
Two major groups of tasks: (1) Tidal analysis, (2) Motion analysis (e.g. Bifrost GPS)
(1) Programs that make use of urtap tidal results:
m/prl2blq.f -> ~/Oload/p/mpp/grdol.f -> visualize,
compare with models
m/urtpp.f -> predict
m/grdcvem.f -> visualize confidence
tsplot, ~/sas/... -> time-series plots and analysis
(2) There is a help file HOW.TO
GPS in which urtap applications appear.
o
Instruction file structure:
===================
1 Open file block
2 Namelist
3 Misc
4 END OF INSTRUCTIONS
codeword
o
Misc consists of individual instruction blocks of
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
END
o
MODEL COVARIANCE
codeword
#n
number of lines that follow
covariance data
...
_______________________________________________________________________________________
o
Defining the "theoretical" signals: Code-word "Tide
Bands" instruction block
====================================================================
Grammar: Underlined characters are functional units
thatocannot be modified.
[ ] indicates
optional code.
____________________________________________________________________________
Order of appearance within "Tide Bands" block:
(1) (...)
TIDE WAVE GROUPS
(...) | <...)
TIDE WAVE GROUPS or INCLUDED TIDE WAVE SUBGROUPS
(2) (~...)
SINGLE-WAVES
o
(3) #...
OBS-DATA
(4) +....
LIN-COMB DATA
____________________________________________________________________________
Select signals for a customized residual:
BEGIN
selectors
END
See the example.
_____________________________________________________*______________________
\ o *
Coding TIDE WAVE GROUPS:
______*______
Identifier = (
(argument selector) = tide symbol
Refer to documentation in WBANDS.f
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
args = (n,tau,s,h,p,N,ps,pi/2)
____________________________________________________________________________
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
____________________________________________________________________________
Coding SINGLE-WAVES:
Identifier = (~
(~ frequency ) = symbol
Frequency in cyc/d
Example:
(~ 4.00d0 ) = S4
____________________________________________________________________________
Coding BOX CARS and SLOPES:
Identifiers = = and %
Box car, i.e. a time-limited offset that is to be estimated:
= ( {From date}{ To date}
) [ comment][ = NAME]
= SYMBOL
Box car 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 1/100 sec}
NAME - a short string (4 chars) to identify the
signal in the result protocol
If not given,
the 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.
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
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.
Equivalence cases for TARGET's 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.
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. This event
is likely
to affect all (TurboRogue) receivers.
____________________________________________________________________________
Designating OBS-DATA= data input into signal matrix using Read_FuF through
Proc_f_on_rec ~/sas/p/urtapu13.f
Identifier = #
#iuninx, trg,rec_mrs,fmt,khms,itz,lddf [ E:tsf_editname] ] = symbol [ -> EOF ]
#iuninx, trg,rec_mrs,fmt,khms,itz,lddf [ R:repairoption]] = symbol [ -> EOF ]
#iuninx, trg,rec_mrs,fmt,khms,itz,[-9
AMP=value] = 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 '(...)' 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 - Call TSF_EDIT (~/sas/p/tsfedit.f)
tsfeditname - char*16 - TSF_EDIT target string
to locate the command section.
R:
- codeword - Call REPAIR (~/sas/p/repair.f)
repairoption - char*4 - option (L to
interpolate linear, else set to DCL)
=
- 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.
lddf -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] 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=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
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. We take the data as is,
khms is ignores 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:
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_UEnd 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, yet some of the subfunctions might like to use for commentary purpose (displays, protocol).
## is synonyme with #_ #: #< #. #$
Here's how to input a series from a multichannel file:
##23 # /geo/hgs/geoc/jpl/jpl_geoc.mc
# 23,'BIN]',-9999.9d0,'L:X_VAL]',0,0,0 R:.] = GEOX
0.001
##23 # /geo/hgs/geoc/jpl/jpl_geoc.mc
# 23,'BIN]',-9999.9d0,'L:Y_VAL]',0,0,0 R:.] = GEOY
0.001
##23 # /geo/hgs/geoc/jpl/jpl_geoc.mc
# 23,'BIN]',-9999.9d0,'L:Z_VAL]',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.
____________________________________________________________________________
Customized residual.
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 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, 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 = +
+iuninx, trg,
rec_mrs, fmt, khms,
itz, lddf, vadd [ E:tsfeditname]
]
Y <- Y + vadd ScaleFactor Filter*R
Again, lddf declares whether scale factor or filter information follows.
vadd - real*8 (sorry) - +1.0d0 | -1.0d0 - add
or subtract the data
from the observations
Other parameters c.f. OBS-DATA BANDS
____________________________________________________________________________
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 bank. Use unit 42 for read-only.
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:
============
Open status is checked. If open,...
41 - Prediction
Error Filter bank of last iteration
12 - eigenvalues,
-vectors and solution parameters:
*.evs-file
13 - tide
band information and response parameters:
*.trs-file
IUN_RESULTS - RSLIST result sheet:
*.prl-file
IUN_RESULTS is a namelist parameter [default=14].
IUMCR - Abridged results
of multi-component analysis output.
IUMCR is a namelist parameter [default=16].
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.
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
.pa.ts predicted
pn + tn
.ra.ts residual
pn + tn
.fa.ts filtered
residual
.sr.ts smoothed
prediction
.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
File name compilation involves replacement of environment parameters.
E.g. the following can be specified in the namelist:
TSFILE_NAMES = '${site}${tapcomp} '
____________________________________________________________________________
PARAMETERS
Array size
NYD - max. time series length.
NBANDS - number of signal components
(2 for every tidal wave group).
MXSPEF - maximum order of PEFilters.
____________________________________________________________________________
NAMELIST &PARAM:
================
Parameter
Meaning
[ default ]
-----------------------------------------------------------------------------
EXEC CONTROL
QBATCH - logical - Avoid all prompts
and graphic. 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 wake up: Q_window_residual
use_eigenvalues - char*16 -
"Select No.of eigenvalues"
e.g.'W 6.0 '
['A']
use_solution_options - char*16 - "Solution options"
[' ']
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]
TIDES
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
['TGP']
Demodulation_Tide - char*32 - Darwin symbol
or argument numbers for a
demodulating frequency. E.g. '2,1,1,-1,0,0,0,0' [' ']
Eff_Average_Time - real*8 -
Assume the average tide during the interval
specified is the effective measurement.
[0.0d0]
Iun_FCP - integer - 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). [T]
Q_always_DC - logical - .true. - Always analyse a DC signal.
[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 to the detected DC-level of the input
signal [0]
SITE AND TIME
XSITE,YSITE - real*8 - Site coordinates, geographic, E-long.,
N-lat. [-]
Idate(3) - integer - Requested date of epoch
for start. A negative
Idate(1) uses the first date of the input file
[3*-9]
Start_h - real*8 - Requested hour
after epoch for first sample [0.d0]
Q_Cut_at_Epoch - log. - Cut input at the Idate and Start_h
If false, this is simply a reference date [.false.]
Ldate(3) - integer - Requested ending date, -9 for all [3*-9]
DT - real*8
- Force sampling interval DT (e.g. to undersample)
Zero value: DT is determined from input data [0.d0]
DTUT - real*8 -
D(TDT-UT) in [s]; BHI-Bull.B: D(TAI-UT) + 32.7 [56.0]
IALTG - integer - Global reference
for grav. and site radius. [2]
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 - Max.duration for
reading procedure.
[(all)]
Ndata_Use - integer - Truncate input data, -1: don't
[-1]
NX_Extra - integer - Additional space for reading
#-bands
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]$']
Q_TSF_Edit - logical - Call TSF_EDIT to work on the input
time series [F]
TSF_Edit_Name-char*32 - name tag of the TSF EDIT block
[]
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']
EIGENVALUE ANALYSIS
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()
especially when the user does not yet know the
eigenvalues/-vectors
[.false.]
q_auto_select_eu - 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 - 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. Try auto_select_eu_margin=0.1
[0.0]
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]
QOUTS_I_SEL - logical - Invert the meaning of the hand selected
(customized) resdiual
[.false.]
MISC.
Q_Ocean - logical - set oceanographic phase
convention
[F]
IALTG - integer - Gravity value
at the site - comp.opt.
[2]
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 ]
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]
DATA NOISE (PRE-)WHITENING
DFF(0:L_Dff)- real*8 - Difference filter, L_Dff = 0 | 1,
used for
PRE-whitening (!)
[1.0,-1.0]
Detrend - char*1 - Detrend the data
('Y'|'N'|' ')
' ': detrend if DFF has nonzero response at freq=0
'N': do not detrend although DFF...
[' ']
LPef - integer - Use
a filter (on file unit 41) 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']
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 [0]
(The idea with m_Burg<=0 is that we request to
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]
Q_Downw - logical - Process a Downweight
block in *.ins file
[F]
NRJD - integer - Reject
diurnals & overtones (S1,S2,...)
[-1]
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
MxMiss_Psp - integer - allow indicated number of missing samples
[0]
Wt_Psp - char*4 - Window type
(cf \SAS\window.f)
['KAIS']
Ls_Psp - integer - Autocovariance
length
[256]
Wdp_Psp - real*8 - Design parameter
for Kaiser-window in resid.psp [1.8]
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-internal defaults imply '+L-U'
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
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.
subroutine proc_downw (Target,Display,b,y,m,n,mu,effw,t0,dt)
=============================================================
/sas/p/downws.f
might be newer. Unmaintained commentary follows:
Processes a "down-weight block" on *.INS file (unit 4)
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.
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.
wrk(mu) - real*8 - work space
____________________________________________________________________________
Description of Down-Weight Block
Line# Instruction
Explanation
-------------------------------------------------------------------------
1 Down-weight <name> <; comment>
Code from column 1, spell exactly this way
(or DOWNWEIG or downweig or Downweig)
name - Optional instruction block identifier
A named block will be processed iff
Target is found in name.
name=* to suppress name search.
<2.1> <RDC > F=oo: lu,trg,recmrs,fmt,khms,itz,lddf
lu=file unit for weight data input; additional
reading and filtering options c.f. urtapu13.f
(sas/p/urtapu13.f)
This line is optional.
RDC - remove DC-level in residual series
before downweighting.
Is often necessary !
oo - operator
/ - divide (else: multiply)
Q - take square root of
file value
Further records may be indicted by filter
options.
2.2 n <RDC ><..W=type <..P=window params>
<..Normal..>> <..T=length>
n - number of records following.
W= - window code word.
type - window type-code (e.g. HANQ).
P= - parameter code word; design
parameters for KAIS DOCH and
delayed taper windows.
Normal - code word. Designates that
window is applied in Normal
Equations only. Else: window
is applied also on residuals.
T= - taper code word
length - cosine taper of down-weighted
section, transition length.
Is within the limits indicated
in the records following.
.. - any text. Total length of
instruction: 80 chars.
3.. yy mm dd hh mm ss l w
n weight records
yy..ss - date, time.
l - length (number of samples).
w - weight value.
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
0 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.
______________________________________________________________________________
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.
Weighting:
==========
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:
------------------------
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
Weights from file ("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,-9,0.
0 Normal
Downweight *
F=/: 21,'TOP] Rwd Y:92',-99999.,'(f12.6,t65,f10.0)',0,0,-9,0.
0 Normal
The first instruction is used in the RMS-residual stage, the second
in the other cases, which are "NormalEq" "Test" 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 Chi^2.
Additional options:
Q - quiet
!DT - impose the time interval from the measurement time series.
RDC - remove DC-level in measurement series before weighting operation.
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 n
external 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.