Thanks for using Compiler Explorer
Sponsors
Jakt
C++
Ada
Analysis
Android Java
Android Kotlin
Assembly
C
C3
Carbon
C++ (Circle)
CIRCT
Clean
CMake
CMakeScript
COBOL
C++ for OpenCL
MLIR
Cppx
Cppx-Blue
Cppx-Gold
Cpp2-cppfront
Crystal
C#
CUDA C++
D
Dart
Elixir
Erlang
Fortran
F#
Go
Haskell
HLSL
Hook
Hylo
ispc
Java
Julia
Kotlin
LLVM IR
LLVM MIR
Modula-2
Nim
Objective-C
Objective-C++
OCaml
OpenCL C
Pascal
Pony
Python
Racket
Ruby
Rust
Snowball
Scala
Solidity
Spice
Swift
LLVM TableGen
Toit
TypeScript Native
V
Vala
Visual Basic
Zig
Javascript
GIMPLE
fortran source #1
Output
Compile to binary object
Link to binary
Execute the code
Intel asm syntax
Demangle identifiers
Verbose demangling
Filters
Unused labels
Library functions
Directives
Comments
Horizontal whitespace
Debug intrinsics
Compiler
AARCH64 gfortran 10.5.0
AARCH64 gfortran 11.4.0
AARCH64 gfortran 12.1.0
AARCH64 gfortran 12.2.0
AARCH64 gfortran 12.3.0
AARCH64 gfortran 13.1.0
AARCH64 gfortran 13.2.0
AARCH64 gfortran 6.4
AARCH64 gfortran 7.3
AARCH64 gfortran 8.2
ARM (32bit) gfortran 10.5.0
ARM (32bit) gfortran 11.4.0
ARM (32bit) gfortran 12.1.0
ARM (32bit) gfortran 12.2.0
ARM (32bit) gfortran 12.3.0
ARM (32bit) gfortran 13.1.0
ARM (32bit) gfortran 13.2.0
ARM (32bit) gfortran 6.4
ARM (32bit) gfortran 7.3
ARM (32bit) gfortran 8.2
LOONGARCH64 gfortran 12.2.0
LOONGARCH64 gfortran 12.3.0
LOONGARCH64 gfortran 13.1.0
LOONGARCH64 gfortran 13.2.0
MIPS gfortran 12.1.0
MIPS gfortran 12.2.0
MIPS gfortran 12.3.0
MIPS gfortran 13.1.0
MIPS gfortran 13.2.0
MIPS gfortran 4.9.4
MIPS gfortran 5.5.0
MIPS gfortran 9.5.0
MIPS64 gfortran 12.1.0
MIPS64 gfortran 12.2.0
MIPS64 gfortran 12.3.0
MIPS64 gfortran 13.1.0
MIPS64 gfortran 13.2.0
MIPS64 gfortran 4.9.4
MIPS64 gfortran 5.5.0
MIPS64 gfortran 9.5.0
MIPS64el gfortran 12.1.0
MIPS64el gfortran 12.2.0
MIPS64el gfortran 12.3.0
MIPS64el gfortran 13.1.0
MIPS64el gfortran 13.2.0
MIPS64el gfortran 4.9.4
MIPS64el gfortran 5.5.0
MIPS64el gfortran 9.5.0
MIPSel gfortran 12.1.0
MIPSel gfortran 12.2.0
MIPSel gfortran 12.3.0
MIPSel gfortran 13.1.0
MIPSel gfortran 13.2.0
MIPSel gfortran 4.9.4
MIPSel gfortran 5.5.0
MIPSel gfortran 9.5.0
POWER gfortran 12.1.0
POWER gfortran 12.2.0
POWER gfortran 12.3.0
POWER gfortran 13.1.0
POWER gfortran 13.2.0
POWER64 gfortran 12.1.0
POWER64 gfortran 12.2.0
POWER64 gfortran 12.3.0
POWER64 gfortran 13.1.0
POWER64 gfortran 13.2.0
POWER64 gfortran trunk
POWER64le gfortran 12.1.0
POWER64le gfortran 12.2.0
POWER64le gfortran 12.3.0
POWER64le gfortran 13.1.0
POWER64le gfortran 13.2.0
POWER64le gfortran trunk
RISC-V 32-bits gfortran (trunk)
RISC-V 32-bits gfortran 12.1.0
RISC-V 64-bits gfortran (trunk)
RISC-V 64-bits gfortran 12.1.0
RISCV (32bit) gfortran 11.4.0
RISCV (32bit) gfortran 12.2.0
RISCV (32bit) gfortran 12.3.0
RISCV (32bit) gfortran 13.1.0
RISCV (32bit) gfortran 13.2.0
RISCV64 gfortran 11.4.0
RISCV64 gfortran 12.2.0
RISCV64 gfortran 12.3.0
RISCV64 gfortran 13.1.0
RISCV64 gfortran 13.2.0
SPARC LEON gfortran 12.2.0
SPARC LEON gfortran 12.3.0
SPARC LEON gfortran 13.1.0
SPARC LEON gfortran 13.2.0
SPARC gfortran 12.2.0
SPARC gfortran 12.3.0
SPARC gfortran 13.1.0
SPARC gfortran 13.2.0
SPARC64 gfortran 12.2.0
SPARC64 gfortran 12.3.0
SPARC64 gfortran 13.1.0
SPARC64 gfortran 13.2.0
flang-trunk (flang-new)
power64 AT12.0
power64 AT13.0
power64le AT12.0
power64le AT13.0
s390x gfortran 12.1.0
s390x gfortran 12.2.0
s390x gfortran 12.3.0
s390x gfortran 13.1.0
s390x gfortran 13.2.0
x86-64 gfortran (trunk)
x86-64 gfortran 10.1
x86-64 gfortran 10.2
x86-64 gfortran 10.3
x86-64 gfortran 10.4
x86-64 gfortran 10.5
x86-64 gfortran 11.1
x86-64 gfortran 11.2
x86-64 gfortran 11.3
x86-64 gfortran 11.4
x86-64 gfortran 12.1
x86-64 gfortran 12.2
x86-64 gfortran 12.3
x86-64 gfortran 13.1
x86-64 gfortran 13.2
x86-64 gfortran 4.9.4
x86-64 gfortran 5.5
x86-64 gfortran 6.3
x86-64 gfortran 7.1
x86-64 gfortran 7.2
x86-64 gfortran 7.3
x86-64 gfortran 8.1
x86-64 gfortran 8.2
x86-64 gfortran 8.3
x86-64 gfortran 8.4
x86-64 gfortran 8.5
x86-64 gfortran 9.1
x86-64 gfortran 9.2
x86-64 gfortran 9.3
x86-64 gfortran 9.4
x86-64 ifort 19.0.0
x86-64 ifort 2021.1.2
x86-64 ifort 2021.10.0
x86-64 ifort 2021.11.0
x86-64 ifort 2021.2.0
x86-64 ifort 2021.3.0
x86-64 ifort 2021.4.0
x86-64 ifort 2021.5.0
x86-64 ifort 2021.6.0
x86-64 ifort 2021.7.0
x86-64 ifort 2021.7.1
x86-64 ifort 2021.8.0
x86-64 ifort 2021.9.0
x86-64 ifx (latest)
x86-64 ifx 2021.1.2
x86-64 ifx 2021.2.0
x86-64 ifx 2021.3.0
x86-64 ifx 2021.4.0
x86-64 ifx 2022.0.0
x86-64 ifx 2022.1.0
x86-64 ifx 2022.2.0
x86-64 ifx 2022.2.1
x86-64 ifx 2023.0.0
x86-64 ifx 2023.1.0
x86-64 ifx 2023.2.1
x86-64 ifx 2024.0.0
Options
Source code
PROGRAM ETERNA CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Program ETERNA, version 3.0 930801 FORTRAN 77. C C C C ===== Version for MS-DOS, LAHEY-compiler. ==================== C C C C Analysis of Earth tide observations. This file has 4677 records. C C C C This program may be compiled (after apropriate modification of C C routine GEOEXT): C C C C under operation system UNIX with SUN FORTRAN-compiler. C C C C under operation system MS-DOS with LAHEY F77L3-compiler version C C 5.10, options /4 /Z1. C C C C under operation system MS-DOS with MICROSOFT FORTRAN 5.0 C C compiler using compiler options /Gt /G2 /Od /AH. The main C C program and the routines stored in file ETSUBS.FOR have to be C C compiled separately using FL /c option and the object modules C C have to be linked subsequently using FL ETERNA ETSUBS. C C C C The program ETERNA allows the adjustment of tidal parameters and C C meteorological parameters from hourly Earth tide observations. C C As observations can be used: C C C C tidal potential, gravity tides, tilt tides, vertical and C C horizontal tidal displacements, vertical strain tides, C C horizontal strain tides, areal strain, shear strain, volume C C strain and ocean tides. C C C C For gravimeter records, a priori amplitude factors for the C C tidal waves within one wave group are used from the WAHR-DEHANT C C elliptical, uniformly rotating oceanless Earth with inelastic C C mantle, liquid outer core and elastic inner core (PREM elastic C C Earth model with mantle dispersion from ZSCHAU and WANG 1987) in C C order to correct for the Earth response to tidal constituents of C C different degree and order within one wave group. C C C C There may be used highpass/bandpass filtered Earth tide and C C meteorological observations applying one of eight different C C lowpass or bandreject symmetrical numerical zero phase filters C C (parameter KFILT = 1...8) or original Earth tide and meteorolo- C C gical observations (parameter KFILT=0). C C C C In case of applying highpass filtering, the observation blocks C C must exceed the filter length. Observation blocks shorter than C C the filter length are counted in the block table, but not used. C C In case of highpass filtering, the Earth tide and meteorological C C observations are splitted into a low frequency part (below 0.7 C C cycle per day) and a high frequency part (above 0.7 cycle per C C day) by numerical filtering. The meteorological parameters and C C the tidal parameters for daily to quarter daily waves may be C C adjusted from the highpass filtered Earth tide and meteorologi- C C cal hourly observations. C C C C In case of no highpass filtering, the tidal parameters may be C C adjusted for long- to short periodic waves. In this case, C C additional TSCHEBYSCHEFF polynomial bias parameters of different C C degree (to be input per block) may be adjusted for each block in C C order to compensate the drift of the Earth tide sensor. The C C polynomial degree per block should not exceed the block length C C divided by the longest period of the tidal waves, for which C C parameters will be adjusted. You should use at least one bias C C parameter per block in order to compensate for offsets. C C C C For a block with N hourly observations, numbered from J=1...N, C C the drift for observation no. J can be computed from the ad- C C justed TSCHEBYSCHEFF polynomial bias parameters DBIAS by C C C C DTN=(DBLE(J-1)-DBLE(N-1)*0.5D0)/(DBLE(N-1)*0.5D0) C C DAK(1)=1.D0 C C DAK(2)=DTN C C DRIFT=DBIAS(1)*DAK(1)+DBIAS(2)*DAK(2) C C DO 4560 I=3,NBIAS C C DAK(I)=2.D0*DTN*DAK(I-1)-DAK(I-2) C C 4560 DRIFT=DRIFT+DBIAS(I)*DAK(I) C C C C Although the adjustment of TSCHEBYSCHEFF polynomials has been C C tested with simulated data up to degree 100, it is recommended C C to use low degree polynomials only, especially when using the C C HANN window for least squares adjustment. The drift approxi- C C mation should in any case be checked by e.g. plotting the C C residuals after adjustment. C C C C The error estimation for the adjusted tidal parameters is done C C in two different ways: C C C C - by the least squares adjustment procedure, neglecting a the C C autocorrelation of the noise, C C C C - by FOURIER-amplitude spectrum of the residuals. C C C C The error estimation by FOURIER.amplitude spectrum of the C C residuals is known to be more realistic than the error estima- C C tion from the least squares adjustment procedure because the C C neglected autocorrelation of the noise within the least squares C C adjustment procedure. C C C C For the least squares adjustment of tidal parameters, there may C C be applied the unity window of the HANN-window for each block C C (e.g. SCHUELLER 1976). In the case of no highpass filtering for C C the adjustment of longperiodic tidal parameters, the application C C of the HANN-window may produce large errors of the drift approxi-C C mation at the start and end of a block. This should be checked C C by plotting the residuals, and in case of problems the appli- C C cation of the unity window may provide a more stable drift C C approximation. C C C C C C Program options: C C ---------------- C C C C There may be used C C - highpass filtering of the observations or drift modelling by C C TSCHEBYSCHEFF polynomials, C C - eight numerical FIR filters of different length and quality, C C - up to five additional observed meteorological parameters, C C - four different tidal potential developments (DOODSON 1921, C C CARTWRIGHT-TAYLER-EDDEN 1973, TAMURA 1987 and BUELLESFELD C C 1985) mya be used, C C - for the least squares adjustment, the unity window or the HANN C C window can be applied for the weights. C C C C Program restrictions: C C --------------------- C C C C The program is restricted to the processing of hourly Earth C C tide readings. The number of hourly readings within one block C C (i.e. without interruption) is not restricted, and up to C C 300 blocks may be processed with the current program version. C C C C The number of wave groups to be analyzed is restricted to 85, C C The number of additional meteorological parameters, for which C C linear regression paraneters may be adjusted, is restricted to C C 5. The total number of unknowns is restricted to 175. C C C C Disc file description: C C ---------------------- C C C C ETERNA.INP: formatted unit, on which the tidal observations C C have to be stored before the execution of program C C ETERNA. C C ETERNA.PRN: unit of formatted printer file. C C ETERNA.RES: unit of formatted residual file. The residuals C C stored on this file may be plotted using programs C C ETERNA.FAR: unit of formatted FOURIER amplitude spectrum file. C C This file may be used by plotprogram PLOTSPEC.FOR. C C PLOTDATA.FOR and PLOTHISTO.FOR. C C ETCPOT.DAT: formatted unit, on which the tidal potential C C development has to be stored before the execution C C of program ETERNA. C C ETCPOT.UFT: unformatted unit, on which the tidal potential C C development will be stored (if it does not exist) C C or from which the tidal potential development will C C be read in order to speed up the computation. C C TAPE12... unformatted scratch unit (residuals). C C TAPE20: unformatted direct access unit for storage of the C C observations. C C C C Used routines: C C -------------- C C C C CHOLIN: computes inverse of normal equation matrix. C C ETASTE: computes astronomical elements. C C ETBUFF: stores data for NC channels in buffer. C C ETFILT: parallel filtering for NC channels stored in buffer. C C ETGCOF: computes geodetic coefficients. C C ETGREI: computes GREGORIAN date. C C ETINPD: reads observations. C C ETJULD: computes JULIAN date. C C ETLOVE: computes latitude dependent elastic parameters. C C ETMUTC: computes difference ET minus UTC or TDT minus UTC. C C ETNLPF: does numerical filter definition. C C ETPOTA: computes tidal waves. C C ETSDER: searches for data errors. C C GEOEXT: computes jobtime. C C C C Loop index description within main program: C C ------------------------------------------- C C C C IF is a loop index running over the filter length (1...NFI).C C IG is a loop index running over the wavegroups (1...NGR).C C IM is a loop index running over the meteoro. param. (1...NF). C C IO is a loop index running over the observations (1...NO). C C JO is a loop index running over the observations (1...NO). C C IU is a loop index running over the unknowns (1...NU). C C JU is a loop index running over the unknbowns (1...NU). C C IW is a loop index running over the tidal waves (1...NW). C C IFR is a loop index running over FOURIER frequencies (1...NFR). C C C C Numerical accuracy: C C ------------------- C C C C The program has been tested on CDC CYBER 990 of RRZN Hannover C C with 15 digits in single precision, on IBM-AT with 15 digits C C in DOUBLE PRECISION, and on a SUN SPARC2 under UNIX with 15 C C digits in DOUBLE PRECISION, and achieved the same results. C C C C Execution time: C C --------------- C C C C The CPU execution time of ETERNA depends mainly on the number of C C Earth tide observations to be processed, the tidal potential to C C be used and the number of tidal parameters to be adjusted. C C The execution time has been measured with three different data C C sets on a number of different processors using under operation C C system UNIX the SUN-FORTRAN compiler, under operation system C C MS-DOS the Microsoft 5.0 FORTRAN compiler (abbreviated to C C MS-FTN5) and the LAHEY F77L3 compiler (abbreviated to LAHEY): C C C C C C operation system: MS-DOS MS-DOS MS-DOS MS-DOS C C compiler: MS-FTN5 MS-FTN5 MS-FTN5 LAHEY5 C C processor: 286/287 386DX/387 386DX/387 386DX/387 C C speed: 12 Mc 16 Mc 20 Mc 20 Mc C C C C sample file days C C C C HAL29901.DAT 63.5 (889)s 540.8 s 429.63 s 340.21 s C C BFL24903.DAT 121.0 (2872)s 1629.2 s 1303.66 s 1022.54 s C C BHTT4002.DAT 122.0 (3741)s 2128.2 s 1695.61 s 1018.21 s C C C C C C operation system: MS-DOS MS-DOS MS-DOS MS-DOS C C compiler: MS-FTN5 LAHEY5 MS-FTN5 LAHEY5 C C processor: 486DX 486DX 486DX2 486DX2 C C speed: 33 Mc 33 Mc 66 Mc 66 Mc C C C C sample file days C C C C HAL29901.DAT 63.5 90.24 s 63.66 s 48.17 s 34.91 s C C BFL24903.DAT 121.0 274.68 s 188.72 s 142.48 s 101.55 s C C BHTT4002.DAT 122.0 356.80 s 234.64 s 185.43 s 124.73 s C C C C C C operation system: UNIX UNIX C C compiler: SUN-F77 SUN-F77 C C processor: SUN SPARC2 SUN SPARC2 C C speed: 40 Mc 40 Mc C C 3 users 1 user C C C C sample file days C C C C HAL29901.DAT 63.5 61.66 s 40.69 s C C BFL24903.DAT 121.0 142.54 s 117.62 s C C BHTT4002.DAT 122.0 223.81 s 154.46 s C C C C C C References : C C ------------ C C C C BUELLESFELD, F.-J. 1985: Ein Beitrag zur harmonischen Darstell- C C ung des gezeitenerzeugenden Potentials. Deutsche Geodaeti- C C sche Kommission, Reihe C, Heft Nr. 314, Muenchen 1985. C C C C CARTWRIGHT, D.E. and R.J. TAYLER 1971: New Computations of Tide C C Generating Potential. The Geophysical Journal, vol. 23 C C no. 1, Oxford 1971. C C C C CARTWRIGHT, D.E. and A.C. EDDEN 1973: Corrected Tables of Tidal C C Harmonics. The Geophysical Journal, vol. 33, no. 3, C C Oxford 1973. C C C C CHOJNICKI, T. 1973: Ein Verfahren zur Erdgezeitenanalyse in An- C C lehnung an das Prinzip der kleinsten Quadrate. Mitteilungen C C aus dem Institut fuer Theoretische Geodaesie der Universi- C C taet Bonn Nr. 15, Bonn 1973. C C C C DEHANT, V. 1987: Tidal Parameters for an Inelastic Earth. C C Physics of the Earth and Planetary Interiors, 49, 97-116, C C 1987. C C C C DOODSON, A.T. 1921: The Harmonic Development of the Tide Gene- C C rating Potential. Proceedings of the Royal Society, Series C C A 100, 306-328, London 1921. Reprint in International C C Hydrographic Revue vol.31 no. 1, Monaco 1954. C C C C PERTSEV, B. 1957: On the calculation of drift curve in observa- C C tion of bodily tides. Bulletin d' Informations, Marees C C Terrestres, no. 5, 71-72, Bruxelles 1957. C C C C PERTSEV, B. 1959: Ob outchetie spolzaniya nulia pir nabloudenij C C ouprougikh prilivov, Izv. Akad. Naouk SSR, no. 4, 1959. C C C C SCHUELLER, K. 1976: Ein Beitrag zur Auswertung von Erdgezeiten- C C registrierungen. Deutsche Geodaetische Kommission, Reihe C, C C no. 227, Muenchen 1976. C C C C TAMURA, Y. 1987: A Harmonic Development of the Tide-generating C C Potential. Bulletin d'Informations Marees Terrestres no. 99,C C 6813-6855, Bruxelles 1987. C C C C WAHR, J.M. 1981: Body tides on an elliptical, rotating, elastic C C and oceanless earth. Geophysical Journal of the Royal C C astronomical Society, vol. 64, 677-703, 1981. C C C C WENZEL, H.-G. 1974: The correction of tidal force development C C to ellipsoidal normal. Bulletin d'Informations Marees C C Terrestres, Vol. 68, 3748-3790, Bruxelles 1974. C C C C WENZEL, H.-G. 1976: Some remarks to the analysismethod of C C Chojnicki. Bulletin d'Informations Marees Terrestres, C C Vol. 73, 4187-4191, Bruxelles 1976. C C C C WENZEL, H.-G. 1976: Zur Genauigkeit von gravimetrischen Erdge- C C zeitenbeobachtunngen. Wissenschaftliche Arbeiten der Lehr- C C stuehle fuer Geodaesie, Photogrammetrie und Kartographie C C an der Technischen Universitaet Hannover NR. 67, Hannover C C 1976. C C C c WENZEL, H.-G. 1977: Estimation of accuracy for the Earth tide C C analysis results. Bulletin d'Informations, Marees C C Terrestres, Vol. 76, 4427-4445, Bruxelles 1977. C C C C WENZEL, H.-G. 1993: Tidal data processing on a PC. Presented C C at 12th International Symposium on Earth Tides, August 4 C C to 8, Beijing 1993. C C C C WILHELM, H. and W. ZUERN 1984: Tidal forcing field. C C In: LANDOLT-BOERNSTEIN, Zahlenwerte und Funktionen aus C C Naturwissenschaften und Technik, New series, group V, Vol. C C 2, Geophysics of the Solid Earth, the Moon and the Planets,C C Berlin 1984. C C C C ZSCHAU, J. and R. WANG 1987: Imperfect elasticity in the Earth's C C mantle. Implications for Earth tides and long period C C deformations. Proceedings of the 9th International Sym- C C posium on Earth Tides, New York 1987, pp. 605-629, editor C C J.T. KUO, Schweizerbartsche Verlagsbuchhandlung, Stuttgart C C 1987. C C C C ZUERN, W. and H. WILHELM 1984: Tides of the solid Earth. C C In: LANDOLT-BOERNSTEIN, Zahlenwerte und Funktionen aus C C Naturwissenschaften und Technik, New series, group V, Vol. C C 2, Geophysics of the Solid Earth, the Moon and the Planets,C C Berlin 1984. C C C C Program creation: 721230 by H.-G. Wenzel, C C Geodaetisches Institut, C C Universitaet Karlsruhe, C C Englerstr. 7, C C D-76128 KARLSRUHE 1, C C Germany. C C Tel: 0049-721-6082307, C C FAX: 0049-721-694552. C C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de C C Last modification: 930801 by H.-G. Wenzel. C C**********************************************************************C IMPLICIT DOUBLE PRECISION (D) IMPLICIT INTEGER*4 (I-N) CHARACTER CTEXT(8)*10,CENDT*10,CUNIT(11)*8,CHEAD(10)*64 CHARACTER CFILE*8,CMODEL(4)*13,CVERS*10,CWIND(2)*5,CBLOCK*10 CHARACTER CPARM*40 DOUBLE PRECISION DCIN(6),DZERO(6) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C The following DIMENSION statement is concerning the number of C C observed meteorological parameters (multi channel input), which C C is restricted to 5 in the current program version (parameter C C MAXNF). C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PARAMETER (MAXNF=5) INTEGER*2 IUNF(MAXNF) CHARACTER CFY1(MAXNF)*8,CFY2(MAXNF)*8 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C The following DIMENSION statements are concerning the number of C C waves of the tidal potential development, which is 1200 in the C C current program version (parameter MAXNW). C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PARAMETER (MAXNW=1214,MAXST=4*MAXNW) INTEGER*4 NRW(MAXNW) DOUBLE PRECISION DTIDST(MAXST),DTHAM(MAXNW),DTHPH(MAXNW), 1 DTHFR(MAXNW),DBODY(MAXNW) EQUIVALENCE (DTIDST(1),DTHAM(1)),(DTIDST(1215),DTHPH(1)), 1 (DTIDST(2429),DTHFR(1)),(DTIDST(3643),DBODY(1)) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C The following DIMENSION statements are concering the number of C C frequencies, at which the FOURIER amplitude spectrum will be C C computed (parameter MAXFR). C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PARAMETER (MAXFR=1300) DOUBLE PRECISION DOC(MAXFR),DOS(MAXFR) EQUIVALENCE (DOC(1),DTIDST(1)),(DOS(1),DTIDST(1301)) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C The following DIMENSION statements are concerning the number of C C wavegroups to be analyzed, which is 85 in the current program C C version (parameter MAXWG). C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PARAMETER (MAXWG=85) INTEGER*2 NA(MAXWG),NE(MAXWG),INA(MAXWG),INE(MAXWG),IUC(MAXWG), 1 IUS(MAXWG) DOUBLE PRECISION DAM(MAXWG),DFR(MAXWG),DGAM(MAXWG),DDPH(MAXWG), 1 DBOD(MAXWG),DFTFD(MAXWG),DFTFP(MAXWG),DMG(MAXWG),DMP(MAXWG), 2 DXA(MAXWG),DYA(MAXWG) CHARACTER CNSY(MAXWG)*4 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C The following DIMENSION statements are concerning the number of C C blocks of data without interruption, which is 300 in the current C C program version (parameter MAXNB). C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PARAMETER (MAXNB=300) INTEGER IRECA(300),IRECE(300),IDATA(300),ITIMA(300),IDATE(300), 1 ITIME(300),IOB(300),NBIAS(300),IUBIAS(300),IFLAG(300) DOUBLE PRECISION DSAPR(300),DSAPO(300),DTLAG(300) CHARACTER CINSTR(300)*10 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C The following dimension statements are concerning the buffer, C C which is used to store data for NBC channels. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DIMENSION DSTOR(6,2596) DIMENSION ITSTOR(2596),IDSTOR(2596),DFL(6),DFH(6) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C The following DIMENSION statement is concering the number of C C unknown parameters to be adjusted, which is 175 in the current C C program version (parameter MAXNU). C C The array DNVEC used for the storage of the normal equation C C system is equivalence to array DSTOR used for the storage of the C C observations during the numerical filtering. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PARAMETER (MAXNU=175,MAXELM=(MAXNU+1)*(MAXNU+2)/2) DOUBLE PRECISION DAK(MAXNU) DOUBLE PRECISION DNVEC(MAXELM),DX(MAXNU) EQUIVALENCE (DNVEC(1),DSTOR(1,1)) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C The following DIMENSION statement is concerning the different C C numerical lowpass filters which can be used. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PARAMETER (MAXFIL=241) DOUBLE PRECISION DLF(MAXFIL) CHARACTER CFILT*12 COMMON /UNITS/ CUNIT,IC2 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C COMMON /CONST/: C C DPI... 3.1415.... C C DPI2... 2.D0*DPI C C DRAD... DPI/180.D0 C C DRO... 180.D0/DPI C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC COMMON /CONST/ DPI,DPI2,DRAD,DRO COMMON /BLOCKR/ IRECA,IRECE,IDATA,ITIMA,IDATE,ITIME,IOB,NBIAS, 1 DSAPR,DSAPO,DTLAG COMMON /BLOCKC/ CINSTR CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C COMMON /STORE/: C C C C DSTOR... array(1...2596,1...6), in which the Earth tide and C C meteorological observations are stored. C C IDSTOR... array(1..2596), in which the date referring to the C C observations is stored. C C ITSTOR... array(1...2596), in which the time referring to C C the observations is stored. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC COMMON /STORE/ DSTOR,IDSTOR,ITSTOR DATA IUN4/4/,IUN5/15/,IUN6/16/,IUN10/10/,IUN11/11/,IUN12/12/, 1 IUN14/14/,IUN20/20/ DATA CVERS/'3.0 930801'/,CWIND/'UNITY','HANN '/ DATA CENDT/'C*********'/ DATA CMODEL/'DOODSON 1921 ','CTED 1973 ','TAMURA 1987 ', 1 'BUELLESF.1985'/ DATA DZERO/6*0.D0/ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Open the files: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC OPEN(UNIT=IUN4, FILE='ETCPOT.DAT',FORM='FORMATTED') OPEN(UNIT=IUN5, FILE='ETERNA.INP',FORM='FORMATTED') OPEN(UNIT=IUN6, FILE='ETERNA.PRN',FORM='FORMATTED') OPEN(UNIT=IUN10,FILE='ETERNA.RES',FORM='FORMATTED') OPEN(UNIT=IUN11,FILE='ETERNA.FAR',FORM='FORMATTED') OPEN(UNIT=IUN12,FORM='UNFORMATTED',STATUS='SCRATCH') OPEN(UNIT=IUN20,ACCESS='DIRECT',STATUS='SCRATCH',RECL=80) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Unit IUN14 will be opened by the call of routine ETPOTA. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC WRITE(IUN6,17000) CVERS WRITE(*,17000) CVERS CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Read data file name CFILE: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC READ(IUN5,17072) CFILE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Read data file header: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1000 READ(IUN5,17070) (CTEXT(I),I=1,8) WRITE(IUN6,17071) (CTEXT(I),I=1,8) WRITE(IUN10,17070) (CTEXT(I),I=1,8) WRITE(IUN11,17070) (CTEXT(I),I=1,8) IF(CTEXT(1).NE.CENDT) GOTO 1000 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DB..ellipsoidal latitude of the station in degree. C C DL..ellipsoidal longitude of the station in degree, positive C C ast of Greenwhich. C C DH..ellipsoidal height of the station in m. C C DG..gravity of the station in m/s**2. C C DAZ... azimuth of the component in degree (only valid for tilt C C and strain). C C C C Parameter IC defines the Earthtide component. C C IC=-1: tidal potential in m**2/s**2. C C IC= 0: vertical tidal acceleration (gravity tide) in nm/s**2 C C (positive downwards). C C IC= 1: horizontal tidal acceleration (tidal tilt) in azimuth C C DAZ in mas = arc sec/1000. C C IC= 2: vertical tidal displacement in mm. C C IC= 3: horizontal tidal displacement in azimuth DAZ in mm. C C IC= 4: vertical tidal strain in 10**-9 = nstr. C C IC= 5: horizontal tidal strain in azimuth DAZ in 10**-9 = nstr. C C IC= 6: areal tidal strain in 10**-9 = nstr. C C IC= 7: shear tidal strain in 10**-9 = nstr. C C IC= 8: volume tidal strain in 10**-9 = nstr. C C IC= 9: ocean tides in millimeter. C C C C Parameter IR is a printout parameter for the tidal developments. C C IR= 0 no printout of the tidal potential development. C C IR =1 printout of the tidal potential development. C C C C IDA IS A PILOTPARAMETER FOR DATA ERROR CHECK. C C IDA =0 NO DATA ERROR CHECK. C C IDA =1 CHECK FOR DATAERRORS. C C C C DATLIM IS THE UPPER LIMIT FOR DATA ERRORS. C C IF THE DATAERRORS EXCEED DATLIM, THE EXECUTION OF THE PROGRAM C C WILL BE TERMINATED. C C C C The initial epoch, defined by ITYI, ITMI, ITDI, ITHI will be C C used to compute the tidal development with amplitudes DTHAM, C C phases DTHPH and frequencies DTHFR at that specific epoch. C C The transformation from UTC to TDT will be carried out at the C C initial epoch. C C C C ITYI... year of initial epoch. C C ITMI... month of initial epoch. C C ITDI... day of initial epoch. C C ITHI... hour of initial epoch. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC READ(IUN5,17004) CPARM,DB WRITE(IUN10,17004) CPARM,DB READ(IUN5,17004) CPARM,DL WRITE(IUN10,17004) CPARM,DL READ(IUN5,17004) CPARM,DH WRITE(IUN10,17004) CPARM,DH READ(IUN5,17004) CPARM,DG WRITE(IUN10,17004) CPARM,DG READ(IUN5,17004) CPARM,DAZ WRITE(IUN10,17004) CPARM,DAZ READ(IUN5,17003) CPARM,IC WRITE(IUN10,17003) CPARM,IC READ(IUN5,17003) CPARM,IR WRITE(IUN10,17003) CPARM,IR READ(IUN5,17003) CPARM,IDA,DATLIM WRITE(IUN10,17003) CPARM,IDA,DATLIM READ(IUN5,17001) CPARM,ITYI,ITMI,ITDI,ITHI WRITE(IUN10,17001) CPARM,ITYI,ITMI,ITDI,ITHI READ(IUN5,17003) CPARM,KFILT WRITE(IUN10,17003) CPARM,KFILT CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C KFILT is a parameter for selecting a specific numerical filter. C C The default for KFILT is 7. C C KFILT= 0: No lowpass filter is used. C C KFILT= 1: lowpass filter of PERTSEV 1957 with 37 coefficients. C C KFILT= 2: lowpass filter of PERTSEV 1959 with 51 coefficients. C C KFILT= 3: lowpass filter with 143 coefficients is used. C C KFILT= 4: lowpass filter with 239 coefficients is used. C C KFILT= 5: lowpass filter with 49 coefficients is used (taken C C from SCHUELLER'S analysis program HYCON-MC). C C KFILT= 6: lowpass filter of 49 hours length (WENZEL 1993). C C KFILT= 7: bandreject filter of 145 hours length (WENZEL 1993). C C KFILT= 8: bandreject filter of 241 hours length (WENZEL 1993). C C C C Parameter IPROBS enables the printout of the original observa- C C tions. C C IPROBS=0: no printout of original observations. C C IPROBS=1: printout of original observations. C C C C Parameter IPRLF enables the printout of lowpass filtered C C observations. Lowpass filtered observations may be used as an C C estimate of the instrumental drift. C C IPRLF= 0 no printout of lowpass filtered observations. C C IPRLF= 1 printout of lowpass filtered observations. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC READ(IUN5,17003) CPARM,IPROBS WRITE(IUN10,17003) CPARM,IPROBS READ(IUN5,17003) CPARM,IPRLF WRITE(IUN10,17003) CPARM,IPRLF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IMODEL = 0: DOODSON 1921 tidal potential will be used. C C IMODEL = 1: CTED 1973 tidal potential will be used. C C IMODEL = 2: TAMURA 1987 tidal potential will be used. C C IMODEL = 3: BUELLESFELD 1985 tidal potential will be used. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC READ(IUN5,17003) CPARM,IMODEL WRITE(IUN10,17003) CPARM,IMODEL CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Read parameter IRIGID for applying rigid Earth amplitude factors.C C IRIGID is set to 1 for strain tides and oceanic tides. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC READ(IUN5,17003) CPARM,IRIGID WRITE(IUN10,17003) CPARM,IRIGID IF(IC.GE.2) IRIGID=1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Read parameter IHANN for applying a HANN-window for the least C C squares adjustment (see ref. SCHUELLER 1976). C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC READ(IUN5,17003) CPARM,IHANN WRITE(IUN10,17003) CPARM,IHANN IF(IHANN.GT.1) IHANN=1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Read parameter IQUICK for stopping the execution after printout C C of the adjusted parameters. In this case, neither the residuals C C nor thier spectrum will be computed. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC READ(IUN5,17003) CPARM,IQUICK WRITE(IUN10,17003) CPARM,IQUICK CALL GEOEXT(IUN6,DEXTIM) CALL ETNLPF(IUN6,KFILT,NFI,DLF,CFILT) NFI2=NFI/2+1 WRITE(IUN6,17062) CFILT,NFI CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute JULIAN date for initial epoch. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DTHI=DBLE(ITHI) CALL ETJULD(IUN6,ITYI,ITMI,ITDI,DTHI,DTUT) DT=(DTUT-2415020.D0)/36525.D0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute the tidal development for the specfic component IC: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IPRINT=1 IF(IR.EQ.1) IPRINT=2 CALL ETPOTA(IUN4,IUN6,IUN14,IPRINT,IMODEL,DB,DL,DH,DG,DAZ,IC, 1 ITYI,ITMI,ITDI,DTHI,DDT0,MAXNW,NRW,DTHAM,DTHPH,DTHFR,DBODY,NW) IC2=IC+2 IPRINT=0 CLOSE(IUN4) CLOSE(IUN14) CALL GEOEXT(IUN6,DEXTIM) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Read comment on 10 records: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC WRITE(IUN6,17011) CVERS WRITE(IUN6,17007) DO 2070 I=1,10 READ(IUN5,17012) CHEAD(I) WRITE(IUN10,17012) CHEAD(I) 2070 WRITE(IUN6,17013) CHEAD(I) WRITE(IUN6,17008) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Write first 3 comment records on unit IUN11: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 2075 I=1,3 2075 WRITE(IUN11,17012) CHEAD(I) WRITE(IUN11,17024) CUNIT(IC2) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Read number of wavegroups NGR and number of meteorological C C parameters NF. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC READ(IUN5,17005) NGR,NF WRITE(IUN10,17005) NGR,NF NC=NF+1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C The number of wavegroups is restricted to MAXWG. C C The number of meteorological parameters is restricted to MAXNF. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(NGR.GT.MAXWG) WRITE(IUN6,17064) MAXWG IF(NGR.GT.MAXWG) WRITE(*,17064) MAXWG IF(NF.GT.MAXNF) WRITE(IUN6,17065) MAXNF IF(NF.GT.MAXNF) WRITE(*,17065) MAXNF IF(NGR.GT.MAXWG.OR.NF.GT.MAXNF) GOTO 15000 JG=1 IF(NGR.EQ.0) GOTO 2085 WRITE(IUN6,17014) CUNIT(IC2) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Read wave groups and frequency transfer function of the instru- C C ment. The wave group numbers are referring to the first column C C of file ETCPOT.DAT, for all tidal potentials which may be used. C C C C Standard wavegroups for Earth tide analysis, depending on the C C recorded time span: C C C C C C group > 1 month > 6 months > 1 year C C from to from to from to C C C C SA - - - - 11 19 C C SSA - - 20 37 20 37 C C MM 38 90 38 90 38 90 C C MF 91 153 91 153 91 153 C C MTM 154 285 154 285 154 285 C C Q1 286 428 286 428 286 428 C C O1 429 488 429 488 429 488 C C M1 489 537 489 537 489 537 C C P1 - - 538 554 538 554 C C S1 - - - - 555 558 C C K1 538 592 555 585 559 576 C C PSI1 - - - - 577 580 C C PHI1 - - - - 581 592 C C J1 593 634 593 634 593 634 C C OO1 635 739 635 739 635 739 C C 2N2 740 839 740 839 740 839 C C N2 840 890 840 890 840 890 C C M2 891 947 891 947 891 947 C C L2 948 987 948 987 948 987 C C S2 988 1121 988 1008 988 1008 C C K2 - - 1009 1121 1009 1121 C C M3 1122 1204 1122 1204 1122 1204 C C M4 1205 1214 1205 1214 1205 1214 C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DMAXPE=1.D30 IF(NGR.EQ.0) GOTO 2086 DMAXPE=0.D0 DO 2080 IG=1,NGR READ(IUN5,17015) NA(JG),NE(JG),CNSY(JG),DFTFD(JG),DFTFP(JG) WRITE(IUN10,17015) NA(JG),NE(JG),CNSY(JG),DFTFD(JG),DFTFP(JG) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Search for rigid Earth tide amplitude DAM, frequency DFR and C C body tide amplitude factor DBOD of the main wave in the group: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DAM(JG)=0.D0 DO 2081 IW=1,NW IF(NRW(IW).LT.NA(JG)) GOTO 2081 INA(JG)=IW GOTO 2082 2081 CONTINUE 2082 CONTINUE INE(JG)=NW DO 2083 IW=1,NW IF(NRW(IW).LE.NE(JG)) GOTO 2083 INE(JG)=IW-1 GOTO 2084 2083 CONTINUE 2084 CONTINUE NAK=INA(JG) NEK=INE(JG) IUC(JG)=2*JG-1 IUS(JG)=2*JG DO 2090 IW=NAK,NEK IF(IRIGID.EQ.1) DBODY(IW)=1.D0 IF(DTHAM(IW).LE.DAM(JG)) GOTO 2100 DAM(JG)=DTHAM(IW) DFR(JG)=DTHFR(IW)*DRO DBOD(JG)=DBODY(IW) 2100 CONTINUE 2090 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Check frequency of main wave in the group. C C If the frequency is less than 10 deg per hour and the filter C C parameter KFILT is not equal to zero, the tidal parameters C C of the group cannot be estimated. Those wave groups are C C automatically eleminated. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(KFILT.NE.0.AND.DFR(JG).LT.10.D0) GOTO 2080 DPER=360.D0/DFR(JG) DMAXPE=DMAX1(DMAXPE,DPER) WRITE(IUN6,17016) JG,NA(JG),NE(JG),DAM(JG),DFR(JG),CNSY(JG), 1 DFTFD(JG),DFTFP(JG),DBOD(JG) JG=JG+1 2080 CONTINUE 2085 CONTINUE NGR=JG-1 2086 WRITE(IUN6,17017) IF(NF.EQ.0) GOTO 2120 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Read names CFY1 and units CFY2 of meteorological data: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 2110 IM=1,NF READ(IUN5,17018) CFY1(IM),CFY2(IM) WRITE(IUN10,17018) CFY1(IM),CFY2(IM) WRITE(IUN6,17019) IM,CFY1(IM),CFY2(IM) 2110 CONTINUE 2120 IF(NF.EQ.0) WRITE(IUN6,17020) CALL GEOEXT(IUN6,DEXTIM) IF(NGR.EQ.0) GOTO 2215 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute frequency gain DFH for numerical highpass filter. C C Array of amplitudes DTHAM contains now body tide amplitude C C DTHAM*DBODY times frequency gain DFH. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC NAK=INA(1) NEK=INE(NGR) DO 2210 IW=NAK,NEK DTM=DLF(NFI2)*0.5D0 DNFI2=DBLE(NFI+1)*0.5D0 DO 2260 J=1,NFI2-1 DF=DBLE(J)-DNFI2 2260 DTM=DTM+DLF(J)*DCOS(DTHFR(IW)*DF) DGAIN=1.D0-2.D0*DTM DTHAM(IW)=DTHAM(IW)*DBODY(IW)*DGAIN 2210 CONTINUE 2215 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Read Earth tide and meteorological observations. C C NC = number of channels. C C NB = number of blocks. C C NREC = number of data records on direct access unit IUN20. C C This program version is restricted to 300 blocks. If you want C C to use more blocks, you have to modify the dimension statements C C associated with parameter MAXNB. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC NC=NF+1 CALL ETINPD(IUN5,IUN6,IUN20,IPROBS,NC,NB,NREC,NERR) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Check all blocks: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 3060 JB=1,NB IF(KFILT.NE.0) NBIAS(JB)=0 IFLAG(JB)=1 IB=IOB(JB) IF(IB.GT.NFI) GOTO 3050 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C The current block is shorter than the filter length and will C C thus be eliminated from the analysis. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC WRITE(IUN6,17086) JB,IDATA(JB),ITIMA(JB),IDATE(JB),ITIME(JB) IFLAG(JB)=0 3050 CONTINUE 3060 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Write observation summary: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 3530 WRITE(IUN6,17011) CVERS WRITE(*,17011) CVERS WRITE(IUN6,17007) WRITE(*,17007) DO 3540 I=1,10 WRITE(IUN6,17013) CHEAD(I) 3540 WRITE(*,17013) CHEAD(I) WRITE(IUN6,17008) WRITE(*,17008) WRITE(IUN6,17029) CUNIT(IC2) WRITE(*,17029) CUNIT(IC2) NOBS=0 DOBSH=0.D0 DO 3550 JB=1,NB NOBS=NOBS+IOB(JB)*IFLAG(NB) DT=DBLE(IOB(JB))/24.0D0 DOBSH=DOBSH+DT IF(KFILT.GT.0) GOTO 3555 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Check number of bias parameters: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(NGR.EQ.0) DMAXPE=2.D0 DMAXNBI=DMAX1(2.D0,DBLE(IOB(JB))/DMAXPE) MAXNBI=DMAXNBI IF(NBIAS(JB).LE.MAXNBI) GOTO 3555 WRITE(IUN6,17076) JB,MAXNBI WRITE(*,17076) JB,MAXNBI 3555 WRITE(IUN6,17030) JB,IDATA(JB),ITIMA(JB),IDATE(JB), 1 ITIME(JB),DT,DSAPR(JB),DTLAG(JB),NBIAS(JB),IFLAG(JB) WRITE(*,17030) JB,IDATA(JB),ITIMA(JB),IDATE(JB), 1 ITIME(JB),DT,DSAPR(JB),DTLAG(JB),NBIAS(JB),IFLAG(JB) 3550 CONTINUE WRITE(IUN6,17031) NB,IDATA(1),ITIMA(1),IDATE(NB), 1 ITIME(NB),DOBSH WRITE(*,17031) NB,IDATA(1),ITIMA(1),IDATE(NB), 1 ITIME(NB),DOBSH DRAY=360.D0/DBLE(NOBS) IF(IHANN.EQ.1) DRAY=2.D0*DRAY DOBSD=DBLE(NOBS)/24.D0 WRITE(IUN6,17032) NOBS*(1+NF) WRITE(*,17032) NOBS*(1+NF) CLOSE(IUN5) CALL GEOEXT(IUN6,DEXTIM) IF(NGR.EQ.0) GOTO 3570 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Check RAYLEIGH-criterion for wavegrouping: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 3560 JG=2,NGR IG=JG-1 DDFR=DFR(JG)-DFR(IG) IF(DABS(DDFR).GE.DRAY) GOTO 3560 WRITE(IUN6,17038) IG,JG WRITE(*,17038) IG,JG 3560 CONTINUE 3570 CONTINUE IF(NERR.EQ.0) GOTO 4000 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Sequence errors ocuured during input of data, program stops: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC WRITE(*,17034) WRITE(IUN6,17034) STOP 4000 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Read data from direct access unit IUN20 and store part of them C C them in buffer (doues routine ETBUFF) : C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC NDL=0 NOB=0 DO 4210 JB=1,NB IF(IFLAG(JB).EQ.0) GOTO 4210 IF(IPRLF.NE.0) WRITE(IUN6,17035) CVERS,JB NO=IOB(JB) IO=0 DO 4220 IREC=IRECA(JB),IRECE(JB) READ(IUN20,REC=IREC) IDAT,ITIM,(DCIN(IC),IC=1,NC) IO=IO+1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Buffering of Earth tide and meteorological observations: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL ETBUFF(NFI,NC,ISTOR,IDAT,ITIM,DCIN,IA,IE) IF(IO.LT.NFI) GOTO 4220 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Lowpass and highpass filtering: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL ETFILT(NFI,NFI2,DLF,NC,IA,IE,IDF,ITF,DFL,DFH) IF(IPRLF.NE.0) WRITE(IUN6,17033) IDF,ITF,(DFL(JC),JC=1,NC) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Store highpass filtered observation vector on IUN20: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IRECF=IREC-NFI2+1 WRITE(IUN20,REC=IRECF) IDF,ITF,(DFH(IC),IC=1,NC) NOB=NOB+1 4220 CONTINUE 4210 CONTINUE C IF(IDA.EQ.0) GOTO 4300 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Search for data errors: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALL ETSDER(IUN6,JB,DAL,MAXOBS,NO,DATLIM,NDLB) C NDL=NDL+NDLB C 4300 CONTINUE C IF(NDL.EQ.0) GOTO 4450 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C LIST TOTAL NUMBER OF DATA ERRORS. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C WRITE(IUN6,17011) CVERS C WRITE(*,17011) CVERS C WRITE(IUN6,17047) NDL,DATLIM,CUNIT(IC2) C WRITE(*,17047) NDL,DATLIM,CUNIT(IC2) C 4450 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute normal equation system sequentially. We use only one C C observation at a time. C C NU = number of unknown parameters. C C The sequence of unknowns is 2*NGR tidal parameters (2 for each C C wave group), NF meteorological regression parameters, and C C bias parameters for each block (in case of no highpass C C filtering). C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC NU=2*NGR+NF DO 4460 IM=1,NF 4460 IUNF(IM)=2*NGR+IM CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Bias unknowns per block for KFILT=0. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(KFILT.GT.0) GOTO 4480 IU=NU+1 DO 4470 JB=1,NB IUBIAS(JB)=IU 4470 IU=IU+NBIAS(JB) NU=IU-1 NDF=NOB-NU 4480 CONTINUE IF(NU.LE.MAXNU) GOTO 4485 WRITE(IUN6,17069) MAXNU,NU WRITE(*,17069) MAXNU,NU 4485 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Initialize normal equation matrix and right hand side: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC NELEM=(NU+1)*(NU+2)/2 DO 4505 K=1,NELEM 4505 DNVEC(K)=0.D0 DO 4600 JB=1,NB IF(IFLAG(JB).EQ.0) GOTO 4600 NO=IOB(JB)-NFI+1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute weight for current block from a priori standard C C deviation DSAPR. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DPB=1.D0/(DSAPR(JB)**2) CALL DATUM(IDATA(JB),ITIMA(JB),ITY,ITM,ITD,ITH) DTH=DBLE(ITH) CALL ETJULD(IUN6,ITY,ITM,ITD,DTH,DT) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Correct for TDT and instrumental time lag: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL ETMUTC(IUN6,IPRINT,DT,DDT) DT=(DT-DTUT)*24.D0+DBLE(NFI-1)/2.D0+(DDT-DDT0-DTLAG(JB))/3600.D0 DIB1=DPI2/DBLE(NO) DIB2=DBLE(NO-1)*0.5D0 DNO2=DBLE(NO-1)*0.5D0 DNO2I=1.D0/DNO2 IREC1=IRECA(JB)+NFI2-1 IREC2=IRECE(JB)-NFI2+1 JO=0 DO 4590 IREC=IREC1,IREC2 READ(IUN20,REC=IREC) IDAT,ITIM,(DFH(JC),JC=1,NC) JO=JO+1 DALJO=DFH(1) DO 4510 IU=1,NU 4510 DAK(IU)=0.D0 IF(NGR.EQ.0) GOTO 4535 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute observation equation and store it in array DAK. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 4530 IG=1,NGR DO 4530 IW=INA(IG),INE(IG) DC=DTHFR(IW)*DT+DTHPH(IW) DAK(IUC(IG))=DAK(IUC(IG))+DTHAM(IW)*DCOS(DC) 4530 DAK(IUS(IG))=DAK(IUS(IG))+DTHAM(IW)*DSIN(DC) 4535 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Observation equation for meteorological regression parameter: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(NF.EQ.0) GOTO 4550 DO 4540 IM=1,NF 4540 DAK(IUNF(IM))=DFH(IM+1) 4550 CONTINUE IF(KFILT.NE.0) GOTO 4570 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Observation equation for polynomial bias parameters. We use C C TSCHEBYSCHEFF polynomials because they are orthogonal with C C respect to normalized time and have excellent numerical C C properties. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(NBIAS(JB).EQ.0) GOTO 4570 DTN=(DBLE(JO-1)-DNO2)*DNO2I IU=IUBIAS(JB) DAK(IU)=1.D0 IF(NBIAS(JB).EQ.1) GOTO 4570 IU=IU+1 DAK(IU)=DTN IF(NBIAS(JB).EQ.1) GOTO 4570 DO 4560 JU=3,NBIAS(JB) IU=IU+1 4560 DAK(IU)=2.D0*DTN*DAK(IU-1)-DAK(IU-2) 4570 DPBJO=DPB CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute HANN-window for the current block: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(IHANN.EQ.1) DPBJO=DPBJO*(1.D0+DCOS((DBLE(JO-1)-DIB2)*DIB1)) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Update normal equation matrix, upper triangle only: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DAK(NU+1)=DALJO K=0 DO 4580 IU=1,NU+1 DO 4580 JU=1,IU K=K+1 4580 DNVEC(K)=DNVEC(K)+DAK(IU)*DAK(JU)*DPBJO 4590 DT=DT+1.D0 4600 CONTINUE CALL GEOEXT(IUN6,DEXTIM) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Solve normal equation system: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL GEOEXT(IUN6,DEXTIM) INV=1 CALL CHOLIN(IUN6,DNVEC,NU,INV,NSING) CALL GEOEXT(IUN6,DEXTIM) IFIRST=NU*(NU+1)/2 DO 5019 IU=1,NU 5019 DX(IU)=DNVEC(IFIRST+IU) DVVP=DNVEC(NELEM) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute standard deviation of weight unit: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC NDF=NOB-NU DMOBP=DSQRT(DVVP/DBLE(NDF)) WRITE(IUN6,17044) DMOBP,NDF WRITE(*,17044) DMOBP,NDF IF(NGR.EQ.0) GOTO 5035 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute adjusted tidal parameters: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 5030 IG=1,NGR DXA(IG)=DX(2*IG-1) DYA(IG)=DX(2*IG) DCC=DSQRT(DXA(IG)**2+DYA(IG)**2) DCS=90.D0 IF(DABS(DXA(IG)).GE.1.D-11)DCS=-DRO*DATAN2(DYA(IG),DXA(IG)) DGAM(IG)=DCC*DBOD(IG)/DFTFD(IG) DDPH(IG)=DCS+DFTFP(IG) 5030 CONTINUE 5035 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Print adjusted parameters: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC WRITE(IUN6,17028) CVERS,CFILE WRITE(*,17028) CVERS,CFILE WRITE(IUN6,17007) WRITE(*,17007) DO 5400 I=1,10 WRITE(IUN6,17013) CHEAD(I) 5400 WRITE(*,17013) CHEAD(I) WRITE(IUN6,17008) WRITE(*,17008) WRITE(IUN6,17039) WRITE(*,17039) NDF=NOB-NU DMOB=DSQRT(DVV/DBLE(NDF)) DMOBP=DSQRT(DVVP/DBLE(NDF)) WRITE(IUN6,17049) (IDATA(JB),ITIMA(JB),IDATE(JB),ITIME(JB), 1 JB=1,NB) WRITE(*,17049) (IDATA(JB),ITIMA(JB),IDATE(JB),ITIME(JB), 1 JB=1,NB) WRITE(IUN6,17040) ITYI,ITMI,ITDI,ITHI WRITE(*,17040) ITYI,ITMI,ITDI,ITHI WRITE(IUN6,17050) DOBSD WRITE(*,17050) DOBSD WRITE(IUN6,17075) CMODEL(IMODEL+1) WRITE(*,17075) CMODEL(IMODEL+1) IF(IC2.EQ.2.AND.IRIGID.EQ.0) WRITE(IUN6,17080) IF(IC2.EQ.2.AND.IRIGID.EQ.0) WRITE(*,17080) IF(IC2.EQ.2.AND.IRIGID.EQ.1) WRITE(IUN6,17081) IF(IC2.EQ.2.AND.IRIGID.EQ.1) WRITE(*,17081) WRITE(IUN6,17082) CWIND(IHANN+1) WRITE(*,17082) CWIND(IHANN+1) WRITE(IUN6,17062) CFILT,NFI WRITE(*,17062) CFILT,NFI IF(NGR.EQ.0) GOTO 5430 WRITE(IUN6,17055) WRITE(*,17055) WRITE(IUN6,17043) CUNIT(IC2) WRITE(*,17043) CUNIT(IC2) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute standard deviation of adjusted tidal parameters: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 5420 IG=1,NGR IU=2*IG-1 JU=2*IG IEL=IU*(IU+1)/2 JEL=JU*(JU+1)/2 DMG(IG)=DMOBP*DSQRT(DXA(IG)**2*DNVEC(IEL)+DYA(IG)**2* 1 DNVEC(JEL)+2.0*DXA(IG)*DYA(IG)*DNVEC(JEL-1))/ 2 DGAM(IG)*DBOD(IG)**2 DMP(IG)=DMOBP*DSQRT(DYA(IG)**2*DNVEC(IEL)+DXA(IG)**2* 1 DNVEC(JEL)+2.0*DXA(IG)*DYA(IG)*DNVEC(JEL-1))*DRO/ 2 DGAM(IG)*DBOD(IG)**2 DAM(IG)=DAM(IG)*DGAM(IG) DSN=DGAM(IG)/DMG(IG) WRITE(IUN6,17046) NA(IG),NE(IG),CNSY(IG),DAM(IG),DSN,DGAM(IG), 1 DMG(IG),DDPH(IG),DMP(IG) WRITE(*,17046) NA(IG),NE(IG),CNSY(IG),DAM(IG),DSN,DGAM(IG), 1 DMG(IG),DDPH(IG),DMP(IG) 5420 CONTINUE 5430 CONTINUE WRITE(IUN6,17044) DMOBP,NDF WRITE(*,17044) DMOBP,NDF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Print adjusted meteorological parameters: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(NF.EQ.0) GOTO 5450 WRITE(IUN6,17058) WRITE(*,17058) DO 5440 IM=1,NF IU=IUNF(IM) IEL=IU*(IU+1)/2 DM=DSQRT(DNVEC(IEL))*DMOBP WRITE(IUN6,17059) IM,DX(IU),DM,CFY1(IM),CUNIT(IC2),CFY2(IM) 5440 WRITE(*,17059) IM,DX(IU),DM,CFY1(IM),CUNIT(IC2),CFY2(IM) 5450 CONTINUE IF(KFILT.NE.0) GOTO 5470 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Print adjusted TSCHEBYSCHEFF polynomial bias parameters: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC WRITE(IUN6,17077) WRITE(*,17077) DO 5460 JB=1,NB IU=IUBIAS(JB) DO 5455 JU=1,NBIAS(JB) IEL=IU*(IU+1)/2 DSBIAS=DMOBP*DSQRT(DNVEC(IEL)) WRITE(IUN6,17078) JB,JU-1,DX(IU),CUNIT(IC2),DSBIAS,CUNIT(IC2) WRITE(*,17078) JB,JU-1,DX(IU),CUNIT(IC2),DSBIAS,CUNIT(IC2) 5455 IU=IU+1 5460 CONTINUE 5470 CONTINUE CALL GEOEXT(IUN6,DEXTIM) IF(IQUICK.EQ.1) GOTO 15000 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute residuals. We compute only one residual at a time. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DVV=0.D0 DVVP=0.D0 NOB=0 DO 6330 JB=1,NB IF(IFLAG(JB).EQ.0) GOTO 6330 DPB=1.D0/(DSAPR(JB)**2) NO=IOB(JB)-NFI+1 NOB=NOB+NO DIB1=DPI2/DBLE(NO) DIB2=DBLE(NO-1)*0.5D0 CALL DATUM(IDATA(JB),ITIMA(JB),ITY,ITM,ITD,ITH) ITH=ITH+NFI2-1 DTH=DBLE(ITH) CALL ETJULD(IUN6,ITY,ITM,ITD,DTH,DT) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Correct for TDT and instrumental time lag DTLAG: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL ETMUTC(IUN6,IPRINT,DT,DDT) DT=(DT-DTUT)*24.D0+(DDT-DDT0-DTLAG(JB))/3600.D0 DSAPO(JB)=0.D0 DMAXR=0.D0 DPBJO=1.D0 WRITE(IUN12) DT DNO2=DBLE(NO-1)/2.D0 DNO2I=1.D0/DNO2 IREC1=IRECA(JB)+NFI2-1 IREC2=IRECE(JB)-NFI2+1 JO=0 CBLOCK='RESIDUALS' WRITE(IUN10,17002) CINSTR(JB),1.D0,DSAPR(JB),DTLAG(JB),NBIAS(JB), 1 CBLOCK WRITE(IUN10,17006) (DZERO(JC),JC=1,NC),0.D0 DO 6150 IREC=IREC1,IREC2 READ(IUN20,REC=IREC) IDAT,ITIM,(DFH(JC),JC=1,NC) JO=JO+1 DALJO=DFH(1) IF(NGR.EQ.0) GOTO 6095 DO 6090 IG=1,NGR DCC=0.D0 DCS=0.D0 DO 6080 IW=INA(IG),INE(IG) DC=DTHFR(IW)*DT+DTHPH(IW) DCC=DCC+DTHAM(IW)*DCOS(DC) 6080 DCS=DCS+DTHAM(IW)*DSIN(DC) DALJO=DALJO-DCC*DX(2*IG-1)-DCS*DX(2*IG) 6090 CONTINUE 6095 CONTINUE IF(NF.EQ.0) GOTO 6110 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Meteorological regression parameters: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 6100 IM=1,NF 6100 DALJO=DALJO-DFH(IM+1)*DX(IUNF(IM)) 6110 CONTINUE IF(KFILT.NE.0) GOTO 6130 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C TSCHEBYSCHEFF polynomial bias parameters: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(NBIAS(JB).EQ.0) GOTO 6130 IU=IUBIAS(JB) DTN=(DBLE(JO-1)-DNO2)*DNO2I DAK(IU)=1.D0 DALJO=DALJO-DAK(IU)*DX(IU) IF(NBIAS(JB).EQ.1) GOTO 6130 IU=IU+1 DAK(IU)=DTN DALJO=DALJO-DAK(IU)*DX(IU) IF(NBIAS(JB).EQ.2) GOTO 6130 DO 6120 JU=3,NBIAS(JB) IU=IU+1 DAK(IU)=2.D0*DTN*DAK(IU-1)-DAK(IU-2) 6120 DALJO=DALJO-DAK(IU)*DX(IU) 6130 CONTINUE DALJO=-DALJO WRITE(IUN12) DALJO CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute HANN-window for least squares adjustment: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(IHANN.EQ.1) DPBJO=(1.D0+DCOS((DBLE(JO-1)-DIB2)*DIB1)) DV2=DALJO**2 DVV=DVV+DPBJO*DV2 DVVP=DVVP+DPB*DPBJO*DV2 DSAPO(JB)=DSAPO(JB)+DPBJO*DV2 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Write highpass filtered data and residuals to IUN10: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC WRITE(IUN10,17033) IDAT,ITIM,(DFH(JC),JC=1,NC),DALJO CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Store maximum residual of this block: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 6140 DABDAL=DABS(DALJO) IF(DABDAL.LT.DMAXR) GOTO 6150 JOMAX=JO DMAXR=DABDAL 6150 DT=DT+1.D0 WRITE(IUN10,17033) 99999999 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Print maximum residual of this block: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL DATUM(IDATA(JB),ITIMA(JB),ITY,ITM,ITD,ITH) ITH=ITH+JOMAX-1+NFI2 CALL ETGREI(ITY,ITM,ITD,ITH) DSAPO(JB)=DSQRT(DSAPO(JB)/DBLE(NO)*DBLE(NOB)/DBLE(NOB-NU)) WRITE(IUN6,17060) JB,ITY,ITM,ITD,ITH,DMAXR,CUNIT(IC2),DSAPO(JB), 1 CUNIT(IC2) 6330 CONTINUE WRITE(IUN10,17033) 88888888 NDF=NOB-NU DMOB=DSQRT(DVV/DBLE(NDF)) DMOBP=DSQRT(DVVP/DBLE(NDF)) WRITE(IUN6,17044) DMOBP,NDF WRITE(*,17044) DMOBP,NDF WRITE(IUN6,17042) DMOB, CUNIT(IC2) WRITE(*,17042) DMOB, CUNIT(IC2) CALL GEOEXT(IUN6,DEXTIM) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C FOURIER-spectrum of residuals, resolution .1 degree per hour. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC REWIND IUN12 NFR=650 IF(NFR.GT.MAXFR) NFR=MAXFR DDFR=0.1D0 DDFRAD=DDFR*DRAD DO 6700 IFR=1,NFR DOC(IFR)=0.D0 6700 DOS(IFR)=0.D0 DO 6750 JB=1,NB IF(IFLAG(JB).EQ.0) GOTO 6750 NO=IOB(JB)-NFI+1 DIB=DBLE(NO) DIB1=DPI2/DIB DIB2=DBLE(NO-1)*0.5D0 READ(IUN12) DT DO 6740 IO=1,NO READ(IUN12) DALJO DALJO=DALJO*(1.D0+DCOS((DBLE(IO-1)-DIB2)*DIB1)) DO 6730 IFR=1,NFR DOMDT=DBLE(IFR-1)*DDFRAD*DT DOC(IFR)=DOC(IFR)+DALJO*DCOS(DOMDT) 6730 DOS(IFR)=DOS(IFR)+DALJO*DSIN(DOMDT) 6740 DT=DT+1.D0 6750 CONTINUE DOB4=2.D0/DBLE(NOB) WRITE(IUN6,17053) CUNIT(IC2) DO 6760 IFR=1,NFR 6760 DOC(IFR)=DSQRT(DOC(IFR)**2+DOS(IFR)**2)*DOB4 DO 6770 IFR=1,NFR,5 IFR11=IFR+4 DOM=DBLE(IFR-1)*DDFR WRITE(IUN11,17052) DOM,(DOC(JFR),JFR=IFR,IFR11) 6770 WRITE(IUN6,17052) DOM,(DOC(JFR),JFR=IFR,IFR11) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Estimation of noise by FOURIER-spectrum of residuals C C DMSE0 is the average from 0.0 to 2.9 deg per hour. C C DMSE1 is the average from 12.0 to 17.9 deg per hour. C C DMSE2 is the average from 26.0 to 31.9 deg per hour. C C DMSE3 is the average from 42.0 to 47.9 deg per hour. C C DMSE4 is the average from 57.0 to 62.9 deg per hour. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DMSE0=0.D0 DO 6775 IFR=1,30 6775 DMSE0=DMSE0+DOC(IFR) DMSE0=DMSE0/30.D0 IF(KFILT.NE.0) DMSE0=9999.9999D0 DMSE1=0.D0 DO 6780 IFR=121,180 6780 DMSE1=DMSE1+DOC(IFR) DMSE1=DMSE1/60.D0 DMSE2=0.D0 DO 6790 IFR=261,320 6790 DMSE2=DMSE2+DOC(IFR) DMSE2=DMSE2/60.D0 DMSE3=0.D0 DO 6800 IFR=421,480 6800 DMSE3=DMSE3+DOC(IFR) DMSE3=DMSE3/60.D0 DMSE4=0.D0 DO 6805 IFR=571,630 6805 DMSE4=DMSE4+DOC(IFR) DMSE4=DMSE4/60.D0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Print adjusted tidal parameters, estimation of noise by FOURIER- C C spectrum of residuals. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC WRITE(IUN6,17028) CVERS,CFILE WRITE(*,17028) CVERS,CFILE WRITE(IUN6,17007) WRITE(*,17007) DO 6810 I=1,10 WRITE(IUN6,17013) CHEAD(I) 6810 WRITE(*,17013) CHEAD(I) WRITE(IUN6,17008) WRITE(*,17008) WRITE(IUN6,17039) WRITE(*,17039) WRITE(IUN6,17049) (IDATA(JB),ITIMA(JB),IDATE(JB),ITIME(JB), 1 JB=1,NB) WRITE(*,17049) (IDATA(JB),ITIMA(JB),IDATE(JB),ITIME(JB), 1 JB=1,NB) WRITE(IUN6,17040) ITYI,ITMI,ITDI,ITHI WRITE(*,17040) ITYI,ITMI,ITDI,ITHI WRITE(IUN6,17050) DOBSD WRITE(*,17050) DOBSD WRITE(IUN6,17075) CMODEL(IMODEL+1) WRITE(*,17075) CMODEL(IMODEL+1) IF(IC2.EQ.2.AND.IRIGID.EQ.0) WRITE(IUN6,17080) IF(IC2.EQ.2.AND.IRIGID.EQ.0) WRITE(*,17080) IF(IC2.EQ.2.AND.IRIGID.EQ.1) WRITE(IUN6,17081) IF(IC2.EQ.2.AND.IRIGID.EQ.1) WRITE(*,17081) WRITE(IUN6,17082) CWIND(IHANN+1) WRITE(*,17082) CWIND(IHANN+1) WRITE(IUN6,17062) CFILT,NFI WRITE(*,17062) CFILT,NFI IF(NGR.EQ.0) GOTO 6840 WRITE(IUN6,17054) DMSE0,CUNIT(IC2),DMSE1,CUNIT(IC2),DMSE2, 1 CUNIT(IC2),DMSE3,CUNIT(IC2),DMSE4,CUNIT(IC2),CUNIT(IC2) WRITE(*,17054) DMSE0,CUNIT(IC2),DMSE1,CUNIT(IC2),DMSE2, 1 CUNIT(IC2),DMSE3,CUNIT(IC2),DMSE4,CUNIT(IC2),CUNIT(IC2) DO 6830 IG=1,NGR DM=DMSE0 IF(DFR(IG).GT.10.D0) DM=DMSE1 IF(DFR(IG).GT.20.D0) DM=DMSE2 IF(DFR(IG).GT.35.D0) DM=DMSE3 IF(DFR(IG).GT.50.D0) DM=DMSE4 DSN=DAM(IG)/DM DMA=DGAM(IG)/DSN DMB=DRO/DSN WRITE(IUN6,17046) NA(IG),NE(IG),CNSY(IG),DAM(IG),DSN,DGAM(IG), 1 DMA,DDPH(IG),DMB WRITE(*,17046) NA(IG),NE(IG),CNSY(IG),DAM(IG),DSN,DGAM(IG), 1 DMA,DDPH(IG),DMB 6830 CONTINUE 6840 CONTINUE WRITE(IUN6,17044) DMOBP,NDF WRITE(IUN6,17042) DMOB, CUNIT(IC2) WRITE(*,17044) DMOBP,NDF WRITE(*,17042) DMOB, CUNIT(IC2) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Print adjusted regression parameters for meteorological obs.: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(NF.EQ.0) GOTO 6860 WRITE(IUN6,17058) WRITE(*,17058) DO 6850 IM=1,NF IU=IUNF(IM) IEL=IU*(IU+1)/2 DM=DSQRT(DNVEC(IEL))*DMOBP WRITE(IUN6,17059) IM,DX(IU),DM,CFY1(IM),CUNIT(IC2),CFY2(IM) 6850 WRITE(*,17059) IM,DX(IU),DM,CFY1(IM),CUNIT(IC2),CFY2(IM) 6860 CONTINUE IF(KFILT.NE.0) GOTO 6880 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Print TSCHEBYSCHEFF polynomial bias parameters: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC WRITE(IUN6,17077) WRITE(*,17077) DO 6870 JB=1,NB IU=IUBIAS(JB) DO 6870 JU=1,NBIAS(JB) IEL=IU*(IU+1)/2 DSBIAS=DMOBP*DSQRT(DNVEC(IEL)) WRITE(IUN6,17078) JB,JU-1,DX(IU),CUNIT(IC2),DSBIAS,CUNIT(IC2) WRITE(*,17078) JB,JU-1,DX(IU),CUNIT(IC2),DSBIAS,CUNIT(IC2) 6870 IU=IU+1 6880 CONTINUE 15000 CONTINUE WRITE(IUN6,17066) WRITE(*,17066) CALL GEOEXT(IUN6,DEXTIM) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Format statements: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 17000 FORMAT( 1' ******************************************************'/ 2' * *'/ 2' * Program ETERNA, Version ',A10,' FORTRAN 77. *'/ 3' * *'/ 4' * Analysis of earthtide observations. *'/ 5' * *'/ 5' * The Geodetic Institute of the University *'/ 6' * Karlsruhe wishes much success when using *'/ 7' * ETERNA. *'/ 6' ******************************************************'//) 17001 FORMAT(A40,5X,4I5) 17002 FORMAT(A10,5X,2F10.4,F10.3,I10,A10) 17003 FORMAT(A40,I10,F10.3) 17004 FORMAT(A40,F10.3) 17005 FORMAT(2I5) 17006 FORMAT(15X,7F10.3) 17007 FORMAT(5X,'#',66('#'),'#') 17008 FORMAT(5X,'#',66('#'),'#') 17011 FORMAT(6X,'Program ETERNA, version ',A10,' FORTRAN 77.'/) 17012 FORMAT(A64) 17013 FORMAT(5X,'# ',A64,' #') 17014 FORMAT(///' Wave groups and parameters for main wave'/ 1' WD delta: a priori WAHR-DEHANT inelastic amplitude factor'// 1' No. from to ampl. frequency ', 2'freq. transf. func. WD delta'/ 3 16X,A8,1X,1X,7H[deg/h],21X,5H[deg]/) 17015 FORMAT(2I5,1X,A4,3X,4F10.4) 17016 FORMAT(3I5,F10.4,F10.6,1X,A6,F9.4,F8.3,F12.5) 17017 FORMAT(///' Meteorological or hydrological parameters :'// 1' no. parameter unit'//) 17018 FORMAT(A8,2X,A8) 17019 FORMAT(I5,2X,A8,2X,A8) 17020 FORMAT(3H NO) 17024 FORMAT(A8) 17028 FORMAT(///6X,'Program ETERNA, version ',A10,' FORTRAN 77,',4X, 1'file: ', A8/) 17029 FORMAT(/' Summary of observation data :'//' Block From', 1 13X,2HTo,12X,' Days Stdv. Time lag Bias Used'/ 2 48X,A8,5X,'SECS'/) 17030 FORMAT(I6,2(2X,I8,1X,I6),F8.2,2F10.3,2I5) 17031 FORMAT(/' Total'//I6,2(2X,I8,1X,I6),F8.2) 17032 FORMAT(/' Total number of observations:',2X,I10/) 17033 FORMAT(I8,1X,I6,7F10.3) 17034 FORMAT(/ 1' **************************************************************'/ 2' ****** Sequence error(s) occured during input of data. *****'/ 3' ****** These errors are listed on print file TAPE6. *****'/ 3' ****** Please check your input data at TAPE5 and retry. *****'/ 4' ****** Program ETERNA finished the execution. *****'/ 5' ****** (It was not successfull, sorry). *****'/ 4' **************************************************************') 17035 FORMAT(' Program ETERNA, version ',A10,' FORTRAN 77.'// 1' Lowpass filtered observations block no.',I4//) 17036 FORMAT(' Program ETERNA, version ',A10,' FORTRAN 77.'// 1' Highpass filtered observations block no.',I4//) 17038 FORMAT(' ***SEPARATION OF GROUP',I4,' AND GROUP',I4, 1' IS DANGEROUS'/) 17039 FORMAT(/6X,'Summary of observation data :'/) 17040 FORMAT(/6x,'Initial epoch for tidal force : ',I4,1H.,I2,1H.,I2, 1 1H.,I2/) 17041 FORMAT(' Program ETERNA, version ',A10,' FORTRAN 77.'// 1' Highpass filtered data and residual (last column)'/ 2' Blok no.: ',I5/) 17042 FORMAT(' Standard deviation :',F10.3,2X,A8) 17043 FORMAT(/' Adjusted tidal parameters :'// 1' from to wave ampl. signal/ ampl.fac. stdv.', 2' phase lead stdv.'/18X,'[',A8,'] noise',19X,' [deg]', 3' [deg]'/) 17044 FORMAT(/' Standard deviation of weight unit:',F10.3/ 1 ' degree of freedom: ',I10) 17046 FORMAT(5X,2I5,1X,A4,F8.3,F8.1,F10.5,F9.5,2F10.4) 17047 FORMAT(//' ***** Total number of',I5,' data errors exceed given'/ 1' ***** limit of ',F10.4,2X,A6,' in all blocks.'//) 17049 FORMAT(6X,I8,I6,3H...,I8,I6,2X,I8,I6,3H...,I8,I6) 17050 FORMAT(6X,'Number of recorded days in total : ',F8.2) 17052 FORMAT(F10.2,5F10.4) 17053 FORMAT(' FOURIER-spectrun of residuals'/' amplitudes in ',A8/ 1' Frequency',7X,3H.00,7X,3H.10,7X,3H.20,7X,3H.30,7X,3H.40/ 2' [deg/hour]'//) 17054 FORMAT(/6X,'Estimation of noise by FOURIER-spectrum of residuals'/ 1 6X,'0.1 cpd band',F10.4,1X,A8,5X,' 1.0 cpd band',F10.4,1X,A8/ 2 6X,'2.0 cpd band',F10.4,1X,A8,5X,' 3.0 cpd band',F10.4,1X,A8/ 3 6X,'4.0 cpd band',F10.4,1X,A8// 4 6X,'adjusted tidal parameters :'// 5 6X,'from to wave ampl. signal/ ampl.fac. stdv.', 6' phase lead stdv.'/19X,A8,' noise',19X,' [deg]', 7' [deg]'/) 17055 FORMAT(/6X,'Estimation of noise by least squares method.'/ 1 6X,'Influence of autocorrelation not considered.'/) 17058 FORMAT(/6X,'Adjusted meteorological or hydrological parameters:'// 1 6X,'no. regr.coeff. stdv. parameter unit'/) 17059 FORMAT(6X,I3,2F12.5,2X,A8,2X,A8,'/',A8) 17060 FORMAT(/' Maximum residual in block no. ',I10/ 1' Date: ',I4,'.',I2,'.',I2,'.',I2,' IS :',F10.3,2X,A8/ 2' RMS residual is :',F10.3,2X,A8/) 17062 FORMAT(6X,'Numerical filter is ',A12,' with ',I3,' coefficients.') 17064 FORMAT(/' ***** You may not use more than',I5,' wavegroups in', 1' this program version.'/ 2' ***** Sorry, you have to modify program ETERNA.'/ 3' ***** Program ETERNA stops the execution.'/) 17065 FORMAT(/ 1' ***** You may not use more than',I5,' meteorological', 1' parameters in this program version.'/ 2' ***** Sorry, you have to modify program ETERNA.'/ 3' ***** Program ETERNA stops the execution.'/) 17066 FORMAT(/ 1' **********************************************'/ 2' * Program ETERNA finished the execution. *'/ 3' * (Hopefully it was successfull). *'/ 4' **********************************************') 17069 FORMAT(/ 1' ***** You may not use more than ',I5,' unknowns in this', 1' program version.'/ 2' ***** The current number of unknowns is ',I5/ 3' ***** Sorry, you have to reduce the number of bias parameters'/ 4' ***** or to modify program ETERNA.'/ 5' ***** Program ETERNA stops the execution.'/) 17070 FORMAT(8A10) 17071 FORMAT(1X,7A10,A9) 17072 FORMAT(10X,A8) 17075 FORMAT(6X,A13,' tidal potential used.') 17076 FORMAT(' *** Number of bias parameters for block ',I5, 1' should not exceed ',I5) 17077 FORMAT(/' Adjusted TSCHEBYSCHEFF polynomial bias', 1' parameters :'// 2 6X,'block degree bias stdv.'//) 17078 FORMAT(2I10,F13.6,1X,A8,F13.6,1X,A8) 17080 FORMAT(6X,'WAHR-DEHANT-ZSCHAU inelastic Earth model used.') 17081 FORMAT(6X,'Rigid Earth model used.') 17082 FORMAT(6X,A5,' window used for least squares adjustment.') 17086 FORMAT(/' Block no. :',I5,' from : ',I8,I6,' to ',I8,I6/ 1' is shorter than the filter length and is thus eliminated.'/) END C BLOCK DATA CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C BLOCK DATA for program ETERNA, version 930703 FTN77. C C C C Routine creation: 930401 by H.-G. Wenzel. C C Geodaetisches Institut, C C Universitaet Karlsruhe, C C Englerstr. 7, C C D-76128 KARLSRUHE 1, C C Germany. C C Tel: 0049-721-6082307, C C FAX: 0049-721-694552. C C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de C C Last modification: 930703 by H.-G. Wenzel. C C**********************************************************************C IMPLICIT DOUBLE PRECISION (D) CHARACTER CUNIT(11)*8 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C COMMON /CONST/: C C DPI... 3.1415.... C C DPI2... 2.D0*DPI C C DRAD... DPI/180.D0 C C DRO... 180.D0/DPI C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC COMMON /CONST/ DPI,DPI2,DRAD,DRO COMMON /UNITS/ CUNIT,IC2 DATA DPI/3.141592653589793D0/,DPI2/6.283185307179586D0/, 1 DRAD/1.745329251994330D-02/,DRO/57.295779513082320D0/ DATA CUNIT/'(m/s)**2','nm/s**2 ',' mas ',' mm ',' mm ', 1' nstr ',' nstr ',' nstr ',' nstr ',' nstr ',' mm '/ END C SUBROUTINE CHOLIN(IUN6,DNV,NU,INV,NSING) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Routine CHOLIN, version 930610 FORTRAN 77. C C C C The routine CHOLIN solves a normal equation system and computes C C the inverse of the normal equation matrix using CHOLESKY's C C method. Routine CHOLIN is an extended version of routine CHOL C C written by R. FORSBERG. C C C C Input parameter description: C C ---------------------------- C C C C IUN6... Unit number of formatted print file. C C DNV... vector, in which the upper triangle of the normal C C equation system to be solved (including the right C C hand side as last column) has to be stored column- C C wise before the execution of routine CHOLIN. C C Right hand side has to be stored starting at C C DNV(NU*(NU+1)/2+1). In the last element of DNV, C C the square sum of observations may be stored. C C After the execution, the solution vector is stored C C as last column of vector DNV, and the square sum C C of residuals is stored as last element of DNV. C C NU... number of unknowns. C C INV... parameter for computing the inverse of the normal C C equation matrix. For INV=0, the unknowns will be C C computed but the inverse of the normal equation C C matrix will not be computed. This option saves C C computation time for those cases, where the inverse C C of the normal equation matrix is not necessary. C C C C Output parameter description: C C ----------------------------- C C C C NSING... number of columns of the normal equation system, C C in which a singulartity has been detected during C C the CHOLESKY factorization. If NSING is greater C C zero, the unknwon vector and the inverse of the C C normal equation matrix may be errorness. C C C C Used routines: C C -------------- C C C C Storage of normal equation system in vector DNV: C C ------------------------------------------------ C C C C The following sketch shows the upper triangle of a normal C C equation system with 4 unknowns stored in vector DNV. C C C C X1 X2 X3 X4 right hand C C side C C C C I DNV( 1) I DNV( 2) I DNV( 4) I DNV( 7) I DNV(11) I C C I DNV( 3) I DNV( 5) I DNV( 8) I DNV(12) I C C I DNV( 6) I DNV( 9) I DNV(13) I C C I DNV(10) I DNV(14) I C C I DNV(15) I C C C C Execution time: C C --------------- C C C C The execution time depends drastically on the number of C C unknowns, it is about proportional to the third power of NU. C C The execution time of routine CHOLIN has been measured on a C C the following processors : C C C C CY990 : Cyber 990 (CDC2) of RRZN Hannover at August 31., 1987. C C C C 486DX2: IBM-AT compatible PC with 66 MHz speed at June 10.,1993. C C MSFOR is Microsoft 5.0 FORTRAN compiler (real mode), C C LAHEY4 is LAHEX F77L3 compiler version 5.10 (1993) C C (protected mode). C C C C SPARC2: SUN SPARC2 with 28.5 MIPS, SUN F77 compiler. C C C C Execution time for solution only (no inversion): C C C C CY990 CY990 CY990 486DX2 486DX2 SPARC 2 C C FTN FTN VFTN 66 MHz 66 Mhz 28.5 MIPS C C NU: OL=LOW OL=HIGH VL=HIGH MSFOR LAHEY5 F77 C C C C 50 0.049 s 0.013 s 0.009 s 0.110 s 0.060 s 0.035 s C C 100 0.343 s 0.079 s 0.040 s 0.770 s 0.220 s 0.090 s C C 200 2.546 s 0.544 s 0.192 s 6.090 s 1.650 s 0.637 s C C 300 8.292 s 1.727 s 0.540 s 20.16 s 5.490 s 2.199 s C C 400 19.559 s 3.076 s 1.048 s - 13.510 s 5.195 s C C 500 37.937 s 7.629 s 1.938 s - 26.910 s 10.066 s C C C C C C Execution time for solution and inversion: C C C C 486DX2 486DX2 SPARC2 C C 66 MHz 66 Mhz 28.5 MIPS C C NU: MSFOR LAHEY4 F77 C C C C 50 0.270 sec 0.110 sec 0.051 s C C 100 2.370 sec 0.550 sec 0.258 s C C 200 18.010 sec 4.780 sec 2.316 s C C 300 60.580 sec 15.760 sec 8.547 s C C 400 40.700 sec 20.820 s C C 500 83.700 sec 42.031 s C C C C Routine creation: 851101 by R. Forsberg, C C Geodetic Institute, C C Gamlehave Allee 22, C C DK-2920 CHARLOTTENLUND, C C Denmark. C C Last modification: 930610 by H.-G. Wenzel, C C Geodaetisches Institut, C C Universitaet Karlsruhe, C C Englerstr. 7, C C D-76128 KARLSRUHE 1, C C Germany. C C Tel.: 0721-6082301. C C FAX: 0721-694552. C C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de C C**********************************************************************C IMPLICIT REAL*8 (A-H,O-Z) DOUBLE PRECISION DNV(1) NSING=0 DO 50 NR=1,NU+1 I=NR*(NR-1)/2 IR=I DO 40 NC=1,NR DSUM=0. IC=NC*(NC-1)/2 I=I+1 NC1=NC-1 DO 30 NP=1,NC1 30 DSUM=DSUM-DNV(IR+NP)*DNV(IC+NP) DCI=DNV(I)+DSUM IF(NR.NE.NC) THEN DNV(I)=DCI/DNV(IC+NC) GOTO 40 ENDIF IF(NR.GT.NU) THEN DNV(I)=DCI GOTO 40 ENDIF IF(DCI.GT.0.D0) THEN DNV(I)=DSQRT(DCI) GOTO 40 ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Singularity in element no. I: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC NSING=NSING+1 DNV(I)=1.D99 WRITE(IUN6,7001) I,NR,NC WRITE(*,7001) I,NR,NC 40 CONTINUE 50 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Back substitution: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 80 NC=NU,1,-1 I=I-1 IR=I IC=NC*(NC+1)/2 DNV(I)=DNV(I)/DNV(IC) DO 70 NP=NC-1,1,-1 IR=IR-1 IC=IC-1 DNV(IR)=DNV(IR)-DNV(I)*DNV(IC) 70 CONTINUE 80 CONTINUE IF(INV.EQ.0) RETURN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute inverse of the normal equation matrix: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IPIV=NU*(NU+1)/2 IND=IPIV DO 100 I=1,NU DIN=1.D0/DNV(IPIV) DNV(IPIV)=DIN MIN=NU KEND=I-1 LANF=NU-KEND IF(KEND) 140,140,110 110 J=IND DO 120 K=1,KEND DSUM=0.D0 MIN=MIN-1 NR=IPIV NC=J DO 130 L=LANF,MIN NC=NC+1 NR=NR+L 130 DSUM=DSUM+DNV(NC)*DNV(NR) DNV(J)=-DSUM*DIN 120 J=J-MIN 140 IPIV=IPIV-MIN 100 IND=IND-1 DO 180 I=1,NU IPIV=IPIV+I J=IPIV DO 180 K=I,NU DSUM=0.D0 NR=J DO 170 L=K,NU NC=NR+K-I DSUM=DSUM+DNV(NR)*DNV(NC) 170 NR=NR+L DNV(J)=DSUM 180 J=J+K RETURN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Format statements: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 7001 FORMAT(/' *****warning from routine CHOLIN, version 930610.'/ 1' *****singularity in element:',I5,' row:',I5,' column:',I5/ 2' *****diagonal element set to 1.D99.'/ 3' *****unknown set to 0.000.'/ 3' *****execution will be continued.'/) END C SUBROUTINE DATUM(IDAT,ITIM,ITY,ITM,ITD,ITH) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Routine DATUM, version 930704 FORTRAN 77. C C C C The routine DATUM converts an 8 digit date IDAT and an 6 digit C C ITIM to year, month, day, hour. C C C C Input parameter description: C C ---------------------------- C C C C IDAT... date in 8 digit form. 19930701 means July 1st, 1993.C C ITIM... time in 6 digit form. 131214 means 13 hours, C C 12 minutes and 14 seconds. C C C C Output parameter description: C C ----------------------------- C C C C ITY... INTEGER year. C C ITM... INTEGER month. C C ITD... INTEGER day. C C ITH... INTEGER hour. C C C C Routine creation: 930701 by H.-G. Wenzel, C C Geodaetisches Institut, C C Universitaet Karlsruhe, C C Englerstr. 7, C C D-76128 KARLSRUHE 1, C C Germany. C C Tel.: 0721-6082301. C C FAX: 0721-694552. C C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de C C Last Modification: 930704 by H.-G.Wenzel. C C**********************************************************************C IDATA=IDAT ITY=IDATA/10000 IDATA=IDATA-ITY*10000 ITM=IDATA/100 IDATA=IDATA-ITM*100 ITD=IDATA ITH=ITIM/10000 RETURN END C SUBROUTINE ETASTE(IUN6,IPRINT,IMODEL,DLON,ITY,ITM,ITD,DTH,DAS, 1 DASP,DDT0) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Routine ETASTE, version 930531 FORTRAN 77. C C C C The routine ETASTE computes the astronomical elements for C C a specific epoch, given in UTC. This routine compiles under C C operation system UNIX using the SUN-FORTRAN compiler and under C C operation system MS-DOS using the MS 5.0 FORTRAN compiler. C C C C Reference: TAMURA, Y. 1987: A Harmonic Development of the Tide- C C generating Potential. Bulletin d'Informations Marees C C Terrestres no. 99, 6813-6855, Bruxelles 1987. C C C C WENZEL, H.-G. 1976: Zur Genauigkeit von gravimetri- C C schen Erdgezeitenbeobachtunngen. Wissenschaftliche C C Arbeiten der Lehrtsuehle fuer Geodaesie, Photogramme- C C trie und Kartographie an der Technischen Universitaet C C Hannover NR. 67, Hannover 1976. C C C C All variables with D as first character are double precision. C C C C Input parameter description: C C ---------------------------- C C C C IUN6... unit number of formatted printout file. C C IPRINT... printout parameter. For IPRINT=0, nothing will C C be printed on unit IUN6. C C IMODEL... parameter for the used tidal potential model. C C IMODEL = 0: DOODSON 1921 tidal potential. C C IMODEL = 1: CTED 1973 tidal potential. C C IMODEL = 2: TAMURA 1987 tidal potential. C C IMODEL = 3: BUELLESFELD 1985 tidal potential. C C DLON... ellipsoidal longitude of the station referring to C C geodetic reference system GRS80 in degree, C C positive east of Greenwhich. C C ITY... integer year of the epoch (e.g. 1988). C C ITM... integer month of the epoch (e.g. January = 1). C C ITD... integer day of the epoch. C C DTH... hour of the epoch in UTC. C C C C Output parameter description: C C ----------------------------- C C C C DAS(1)... mean local Moontime in degree. C C DAS(2)... mean longitude of the Moon in degree. C C DAS(3)... mean longitude of the Sun in degree. C C DAS(4)... mean longitude of the perigee of the Moon's orbit C C in degree. C C DAS(5)... neagtive mean longitude of the ascending node of C C the Moon's orbit in degree. C C DAS(6)... mean longitude of the perigee of the Suns's orbit C C in degree. C C DAS(7)... argument of Jupiter's opposition in degree (for C C TAMURA's 1987 tidal potential development). C C DAS(8)... argument of Venus's conjunction in degree (for C C TAMURA's 1987 tidal potential development). C C DASP(1...8): time derivatives of the corresponding variables C C DAS in degree per hour. C C DDT0... difference TDT minus UTC at initial epoch in sec. C C C C Used routines: C C -------------- C C C C ETEXTI: computes jobtime. C C ETJULD: computes Julian date. C C C C Numerical accuracy: C C ------------------- C C C C The routine has been tested on IBM Pc with 15 digits accuracy C C in double precision using different compilers. C C C C Execution time: C C --------------- C C C C The execution time is about 0.00017 s per call of routine ETASTE C C on ICM Pc 80486 DX2 with 66 MHz speed. C C C C Routine creation: 880130 BY H.-G.WENZEL. C C Geodaetisches Institut, C C Universitaet Karlsruhe, C C Englerstr. 7, C C D-76128 KARLSRUHE 1, C C Germany. C C Tel.: 0721-6082307, C C FAX: 0721-694552. C C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de C C Last modification: 930531 by H.-G. Wenzel. C C**********************************************************************C IMPLICIT DOUBLE PRECISION (D) IMPLICIT INTEGER*4 (I-N) DOUBLE PRECISION DAS(8),DASP(8) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C COMMON /CONST/: C C DPI... 3.1415.... C C DPI2... 2.D0*DPI C C DRAD... DPI/180.D0 C C DRO... 180.D0/DPI C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC COMMON /CONST/ DPI,DPI2,DRAD,DRO CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute Julian date for epoch ITY, ITM, ITD, DTH. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL ETJULD(IUN6,ITY,ITM,ITD,DTH,DTUJD) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute Universal Time epoch DTUT in Julian centuries referring C C to 1. January 1900. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DTUT=(DTUJD-2415020.0D0)/36525.0D0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute Universal Time epoch DTUT20 in Julian Centuries C C referring to 1. January 2000 12 h. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DTUT20=(DTUJD-2451545.0D0)/36525.0D0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute Terrestrial Dynamical Time epoch DT in Julian Centuries C C referring to 1. January 1990 and DTDT20 referring to 1 January C C 2000 12 h. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL ETMUTC(IUN6,IPRINT,DTUJD,DDT0) DT=DTUT+DDT0/3155760000.0D0 DTDT20=DTUT20+DDT0/3155760000.0D0 IF(IPRINT.GT.0) WRITE(IUN6,7001) ITY,ITM,ITD,DTH,DTUJD,DDT0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute astronomical elements for initial epoch from C C NEWCOMB's formula for the Sun and from BROWN's formulas for the C C Moon (see Astronomical Ephemeris, Explanatory Supplement). C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DT2=DT*DT DT3=DT2*DT IF(IMODEL.EQ.2) GOTO 100 DAS(2)=270.434164D0 +481267.8831417D0*DT-0.0011333D0*DT2 1 +0.0000019D0*DT3 DAS(3)=279.696678D0 +36000.768925D0*DT +0.0003025D0*DT2 DAS(4)=334.329556D0 +4069.0340333D0*DT -0.010325D0*DT2 1 -0.0000125D0*DT3 DAS(5)=100.816725D0 +1934.1420083D0*DT -0.0020778D0*DT2 1 -0.0000022D0*DT3 DAS(6)=281.220833D0 +1.719175D0*DT +0.0004528D0*DT2 1 +0.0000033D0*DT3 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute Jupiter's and Venus's arguments from TAMURA's 1987 C C formulas. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DAS(7)=248.1D0+32964.47D0*DTDT20 DAS(8)= 81.5D0+22518.44D0*DTDT20 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute speeds in degree per hour: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DASP(2)=0.54901652195037D0 -2.58575D-9*DT +6.46D-12*DT2 DASP(3)=0.04106863897444D0 +6.902D-10*DT DASP(4)=0.00464183667960D0 -2.355692D-8*DT-4.278D-11*DT2 DASP(5)=0.00220641342494D0 -4.74054D-9*DT -7.60D-12*DT2 DASP(6)=0.00000196118526D0 +1.03303D-9*DT +1.141D-11*DT2 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute Jupiter's and Venus's speed from TAMURA's 1987 C C formulas. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DASP(7)=32964.47D0/(24.D0*36525.D0) DASP(8)=22518.44D0/(24.D0*36525.D0) DASP(1)=DASP(3)-DASP(2)+15.0D0 DO 10 I=2,8 DAS(I)=DMOD(DAS(I),360.0D0) IF(DAS(I).LT.0.D0) DAS(I)=DAS(I)+360.D0 10 CONTINUE DAS(1)=DAS(3)-DAS(2)+DLON+DTH*15.0D0 IF(DAS(1).LT.0.D0) DAS(1)=DAS(1)+360.0D0 GOTO 200 100 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute astronomical elements from TAMURA's 1987 formulas: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DT2=DTUT20*DTUT20 DT3=DT2*DTUT20 DAL=280.4606184D0 + 36000.7700536D0*DTUT20 + 0.00038793D0*DT2 1 -0.0000000258D0*DT3 DALP=(36000.7700536D0 +2.0D0*0.00038793D0*DTUT20 1 -3.0D0*0.0000000258D0*DT2)/(24.0D0*36525.D0) DT2=DTDT20**2 DS=218.316656D0+481267.881342D0*DTDT20-0.001330D0*DT2 DSP=(481267.881342D0-2.0D0*0.001330D0*DTDT20)/(24.D0*36525.0D0) DH=280.466449D0+36000.769822D0*DTDT20+0.0003036D0*DT2 DHP=(36000.769822D0+2.0D0*0.0003036D0*DTDT20)/(24.D0*36525.0D0) DDS=0.0040D0*DCOS((29.D0+133.0D0*DTDT20)*DRAD) DDSP=(-0.0040D0*133.0D0*DRAD*DSIN((29.D0+133.0D0*DTDT20)*DRAD))/ 1 (24.0D0*36525.0D0) DDH=0.0018D0*DCOS((159.D0+19.D0*DTDT20)*DRAD) DDHP=(-0.0018D0*19.0D0*DRAD*DSIN((159.D0+19.D0*DTDT20)*DRAD))/ 1 (24.0D0*36525.0D0) DAS(1)=DAL-DS+DLON+DTH*15.0D0 DAS(2)=DS+DDS DAS(3)=DH+DDH DAS(4)=83.353243D0 + 4069.013711D0*DTDT20 -0.010324D0*DT2 DAS(5)=234.955444D0 +1934.136185D0*DTDT20 -0.002076D0*DT2 DAS(6)=282.937348D0 + 1.719533D0*DTDT20 +0.0004597D0*DT2 DAS(7)=248.1D0+32964.47D0*DTDT20 DAS(8)= 81.5D0+22518.44D0*DTDT20 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute the speedsin degree per hour: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DASP(1)=DALP-DSP+15.0D0 DASP(2)=DSP+DDSP DASP(3)=DHP+DDHP DASP(4)=(4069.013711D0-2.0D0*0.010324D0*DTDT20)/(24.0D0*36525.0D0) DASP(5)=(1934.136185D0-2.0D0*0.002076D0*DTDT20)/(24.0D0*36525.0D0) DASP(6)=(1.719533D0+2.0D0*0.0004597D0*DTDT20)/(24.0D0*36525.0D0) DASP(7)=32964.47D0/(24.D0*36525.D0) DASP(8)=22518.44D0/(24.D0*36525.D0) DO 110 I=1,8 DAS(I)=DMOD(DAS(I),360.0D0) IF(DAS(I).LT.0.D0) DAS(I)=DAS(I)+360.0D0 110 CONTINUE 200 CONTINUE IF(IPRINT.EQ.0) RETURN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Print astronomical elements: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(IMODEL.NE.2) WRITE(IUN6,7003) (DAS(K),DASP(K),K=1,8) IF(IMODEL.EQ.2) WRITE(IUN6,7004) (DAS(K),DASP(K),K=1,8) 5000 CONTINUE WRITE(IUN6,7030) WRITE(*,7030) RETURN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Format statements: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 7001 FORMAT(' Routine ETASTE, version 930531 FTN 77.'// 1' Astronomic elements for initial epoch ',I4,2I3,F3.0,' UTC'/ 2' Julian date : ',F15.4/ 3' Correction to UTC to give TDT : ',F15.3,' s'/) 7003 FORMAT(' BROWNs formulas for the Moon, NEWCOMBs formulas for the', 1' Sun'/' TAMURAs formulas for Jupiter and Venus'// 2' TAU',F20.11,' deg TAU.',F20.11,' deg/hour'/ 3' S ',F20.11,' deg S. ',F20.11,' deg/hour'/ 4' H ',F20.11,' deg H. ',F20.11,' deg/hour'/ 5' P ',F20.11,' deg P. ',F20.11,' deg/hour'/ 6' N ',F20.11,' deg N. ',F20.11,' deg/hour'/ 7' P1 ',F20.11,' deg P1. ',F20.11,' deg/hour'/ 8' JU ',F20.11,' deg JU. ',F20.11,' deg/hour'/ 9' VE ',F20.11,' deg VE. ',F20.11,' deg/hour'/) 7004 FORMAT(' TAMURAs 1987 formulas are used.'// 1' F1 ',F20.11,' deg F1. ',F20.11,' deg/hour'/ 2' F2 ',F20.11,' deg F2. ',F20.11,' deg/hour'/ 3' F3 ',F20.11,' deg F3. ',F20.11,' deg/hour'/ 4' F4 ',F20.11,' deg F4. ',F20.11,' deg/hour'/ 5' F5 ',F20.11,' deg F5. ',F20.11,' deg/hour'/ 6' F6 ',F20.11,' deg F6. ',F20.11,' deg/hour'/ 7' F7 ',F20.11,' deg F7. ',F20.11,' deg/hour'/ 8' F8 ',F20.11,' deg F8. ',F20.11,' deg/hour'/) 7030 FORMAT(///' ***** Routine ETASTE finished the execution.'/) END C SUBROUTINE ETBUFF(NFI,NC,ISTOR,IDAT,ITIM,DCIN,IA,IE) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Routine ETBUFF, version 930629 FORTRAN 77. C C C C The routine ETBUFF stores ISTOR data for NC channels in buffer. C C The maximum number of channels is 6, the maximum number of data C C per channel is 2523. C C C C Input parameter description: C C ---------------------------- C C C C NFI... filter length. NFI is restricted to be less or C C equal to 2523. C C NC... number of channels to be used. NC is restricted to C C be less or equal to 6. C C ITSTOR... position, at which the new observation vector will C C be stored. C C IDAT... date of the new observation, which will be stored C C in IDSTOR at position ISTOR. C C ITIM... time of the new observation vector, which will be C C stored in ITSTOR at position ITSTOR. C C DCIN... new observation vector(1...6), which will be C C stored in array DSTOR at position ISTOR. C C C C Output parameter description: C C ----------------------------- C C C C IA... start index (first position) for the actual NFI C C elements stored in arrays IDSTOR, ITSTOR, DSTOR. C C IE... end index (last position) for the actual NFI C C elements stored in arrays IDSTOR, ITSTOR, DSTOR. C C C C COMMON /STORE/: C C -------------- C C C C DSTOR... array(1...2596,1...6), in which the Earth tide and C C meteorological observations are stored. C C IDSTOR... array(1..2596), in which the date referring to the C C observations is stored. C C ITSTOR... array(1...2596), in which the time referring to C C the observations is stored. C C C C Routine creation: 910921 by H.-G. Wenzel, C C Geodaetisches Institut, C C Universitaet Karlsruhe, C C Englerstr. 7, C C D-76128 KARLSRUHE 1, C C Germany. C C Tel.: 0721-6082307, C C FAX: 0721-694552. C C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de C C Last modification: 930629 by H.-G. Wenzel. C C**********************************************************************C IMPLICIT DOUBLE PRECISION (D) DIMENSION DSTOR(6,2596),DCIN(6) DIMENSION ITSTOR(2596),IDSTOR(2596) COMMON /STORE/ DSTOR,IDSTOR,ITSTOR DATA MAXSTOR/2596/ ISTOR=ISTOR+1 IF(ISTOR.GT.MAXSTOR) GOTO 1000 IDSTOR(ISTOR)=IDAT ITSTOR(ISTOR)=ITIM DO 10 J=1,NC 10 DSTOR(J,ISTOR)=DCIN(J) IE=ISTOR IA=ISTOR-NFI+1 RETURN 1000 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Shift the last NFI elements of arrays IDSTOR, ITSTOR and DSTOR C C to their beginnings. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 1010 I=1,NFI IDSTOR(I)=IDSTOR(MAXSTOR-NFI+I) ITSTOR(I)=ITSTOR(MAXSTOR-NFI+I) DO 1010 J=1,NC 1010 DSTOR(J,I)=DSTOR(J,MAXSTOR-NFI+I) ISTOR=NFI ISTOR=ISTOR+1 IDSTOR(ISTOR)=IDAT ITSTOR(ISTOR)=ITIM DO 1020 J=1,NC 1020 DSTOR(J,ISTOR)=DCIN(J) IE=ISTOR IA=ISTOR-NFI+1 RETURN END C SUBROUTINE ETFILT(NFI,NFI2,DFIL,NC,IA,IE,IDF,ITF,DFL,DFH) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Routine ETFILT, version 930703 FORTRAN 77. C C C C The routine does numerical filtering of the data stored in C C array DSTOR using a symmetrical numerical FIR lowpass filter. C C C C All parameters with D as first character are DOUBLE PRECISION. C C C C Input parameter description: C C ---------------------------- C C C C NFI... number of filter coefficients (filter length). C C NFI2... index of central filter coefficient. C C DFIL... array of lowpass filter coefficients (1...NFI). C C C C Output parameter description: C C ----------------------------- C C C C IDF... INTEGER DATE OF FILTERED DATA. C C ITF... TIME OF FILTERED DATA. C C DFL... lowpass filtered observation vector (1...6). C C DFH... highpass filtered observation vector (1...6). C C C C COMMON /STORE/: C C -------------- C C C C DSTOR... array(1...2596,1...6), in which the Earth tide and C C meteorological observations are stored. C C IDSTOR... array(1..2596), in which the date referring to the C C observations is stored. C C ITSTOR... array(1...2596), in which the time referring to C C the observations is stored. C C C C Used routines: none. C C -------------- C C C C Routine creation: 910629 by H.-G. Wenzel, C C Geodaetisches Institut, C C Universitaet Karlsruhe, C C Englerstr. 7, C C D-76128 KARLSRUHE 1, C C Germany. C C Tel.: 0721-6082307, C C FAX: 0721-694552. C C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de C C Last modification: 930703 by H.-G. Wenzel. C C**********************************************************************C IMPLICIT DOUBLE PRECISION (D) DIMENSION DFIL(241),DSTOR(6,2596),DFL(6),DFH(6) DIMENSION ITSTOR(2596),IDSTOR(2596) COMMON /STORE/ DSTOR,IDSTOR,ITSTOR DO 10 J=1,NC 10 DFL(J)=DSTOR(J,IA+NFI2-1)*DFIL(NFI2) DO 20 I=1,NFI2-1 DO 20 J=1,NC 20 DFL(J)=DFL(J)+(DSTOR(J,IA+I-1)+DSTOR(J,IE+1-I))*DFIL(I) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute highpass filtered observations by subtracting the C C lowpass filtered observation from the original observation: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 30 J=1,NC 30 DFH(J)=DSTOR(J,IA+NFI2-1)-DFL(J) IDF=IDSTOR(IA+NFI2-1) ITF=ITSTOR(IA+NFI2-1) RETURN END C SUBROUTINE ETGCOF(IUN6,IPRINT,DLAT,DLON,DH,DGRAV,DAZ,IC,DGK,DPK) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Routine ETGCOF, version 930706 FORTRAN 77. C C C C The routine ETGCOF computes the geodetic coefficients for C C the tidal potential developments, DOODSON's normalization. C C C C All variables with D as first character are DOUBLE PRECISION. C C C C Input parameter description: C C ---------------------------- C C C C IUN6... formatted line printer unit. C C DLAT... ellipsoidal latitude in degree, referring to C C geodetic reference system GRS80. C C DLON... ellipsoidal longitude in degree, referring to C C geodetic reference system GRS80, positiv east of C C Greenwhich. C C DH... ellipsoidal height in meter, referring to geodetic C C reference system GRS80. C C DGRAV... gravity in m/s**2. If DGRAV less than 9.50 m/s**2, C C DGRAV will be overwritten by normal gravity C C referring to geodetic reference system 1980. C C DAZ... azimuth in degree from north direction counted C C clockwise (necessary for tidal tilt only). C C IC... Earth tide component to be computed. C C IC=-1: tidal potential, geodetic coefficients C C in m**2/s**2. C C IC= 0: vertical tidal acceleration (gravity tide), C C geodetic coefficients in nm/s**2 (positive C C down). C C IC= 1: horizontal tidal acceleration (tidal tilt) C C in azimuth DAZ, geodetic coefficients in C C mas = arc sec/1000. C C IC= 2: vertical tidal displacement, geodetic C C coefficients in mm. C C IC= 3: horizontal tidal displacement in azimuth C C DAZ, geodetic coefficients in mm. C C IC= 4: vertical tidal strain, geodetic coefficients C C in 10**-9 = nstr. C C IC= 5: horizontal tidal strain in azimuth DAZ, C C geodetic coefficients in 10**-9 = nstr. C C IC= 6: areal tidal strain, geodetic coefficients C C in 10**-9 = nstr. C C IC= 7: shear tidal strain, geodetic coefficients C C in 10**-9 = nstr. C C IC= 8: volume tidal strain, geodetic coefficients C C in 10**-9 = nstr. C C IC= 9: ocean tides, geodetic coefficients in C C millimeter. C C C C Output parameter description: C C ----------------------------- C C C C DGK... array (1...12) of geodetic coefficients. C C DGK(1)... geodetic coefficient for degree 2 and order 0. C C DGK(2)... geodetic coefficient for degree 2 and order 1. C C DGK(3)... geodetic coefficient for degree 2 and order 2. C C DGK(4)... geodetic coefficient for degree 3 and order 0. C C DGK(5)... geodetic coefficient for degree 3 and order 1. C C DGK(6)... geodetic coefficient for degree 3 and order 2. C C DGK(7)... geodetic coefficient for degree 3 and order 3. C C DGK(8)... geodetic coefficient for degree 4 and order 0. C C DGK(9)... geodetic coefficient for degree 4 and order 1. C C DGK(10)... geodetic coefficient for degree 4 and order 2. C C DGK(11)... geodetic coefficient for degree 4 and order 3. C C DGK(12)... geodetic coefficient for degree 4 and order 4. C C C C DPK... array (1...12) of phases in degree. C C DPK(1)... phase in degree for degree 2 and order 0. C C DPK(2)... phase in degree for degree 2 and order 1. C C DPK(3)... phase in degree for degree 2 and order 2. C C DPK(4)... phase in degree for degree 3 and order 0. C C DPK(5)... phase in degree for degree 3 and order 1. C C DPK(6)... phase in degree for degree 3 and order 2. C C DPK(7)... phase in degree for degree 3 and order 3. C C DPK(8)... phase in degree for degree 4 and order 0. C C DPK(9)... phase in degree for degree 4 and order 1. C C DPK(10).. phase in degree for degree 4 and order 2. C C DPK(11).. phase in degree for degree 4 and order 3. C C DPK(12).. phase in degree for degree 4 and order 4. C C C C Used routines: C C -------------- C C C C ETLOVE: computes latitude dependent elastic parameters. C C C C Numerical accuracy: C C ------------------- C C C C The routine has been tested under operation system MS-DOS and C C UNIX in double precision (8 byte words = 15 digits) using C C different compilers. C C C C References: C C C C WILHELM, H. and W. ZUERN 1984: Tidal forcing field. C C In: LANDOLT-BOERNSTEIN, Zahlenwerte und Funktionen aus C C Naturwissenschaften und Technik, New series, group V, C C Vol. 2, Geophysics of the Solid Earth, the Moon and the C C Planets, Berlin 1984. C C C C ZUERN, W. and H. WILHELM 1984: Tides of the solid Earth. C C In: LANDOLT-BOERNSTEIN, Zahlenwerte und Funktionen aus C C Naturwissenschaften und Technik, New series, group V, Vol. C C 2, Geophysics of the Solid Earth, the Moon and the Planets,C C Berlin 1984. C C C C Routine creation: 880129 by H.-G. Wenzel, C C Geodaetisches Institut, C C Universitaet Karlsruhe, C C Englerstr. 7, C C D-76128 KARLSRUHE 1, C C Germany. C C Tel.: 0721-6082307, C C FAX: 0721-694552. C C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de C C Last modification: 930706 by H.-G. Wenzel. C C**********************************************************************C IMPLICIT DOUBLE PRECISION (D) CHARACTER CUNIT(11)*8 DIMENSION DGK(12),DPK(12),DGX(12),DGY(12),DGZ(12) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C COMMON /LOVE/ contains gravimeter factors, LOVE-numbers, SHIDA- C C numbers and tilt factors for degree 2...4 at latitude DLAT: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DIMENSION DGLAT(12),DHLAT(12),DKLAT(12),DLLAT(12),DTLAT(12) COMMON /LOVE/ DOM0,DOMR,DGLAT,DGR,DHLAT,DHR,DKLAT,DKR,DLLAT,DLR, 1 DTLAT,DTR CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C COMMON /CONST/: C C DPI... 3.1415.... C C DPI2... 2.D0*DPI C C DRAD... DPI/180.D0 C C DRO... 180.D0/DPI C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC COMMON /CONST/ DPI,DPI2,DRAD,DRO COMMON /UNITS/ CUNIT,IC2 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Definition of parameters of Geodetic Reference System 1980. C C DEA is major semi axis in meter. C C DEE is square of first excentricity (without dimension). C C DEGM is geocentric gravitational constant in m*3/s**2. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DATA DEA/6378137.00D0/,DEE/6.69438002290D-3/,DEGM1/398600.5D0/ DEGM=DEGM1*1.D9 IF(IPRINT.GT.0) WRITE(IUN6,7000) DEA,DEE,DEGM1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DCLAT is cos and DSLAT is sin of ellipsoidal latitude. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DCLAT=DCOS(DLAT*DRAD) DSLAT=DSIN(DLAT*DRAD) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute normal gravity in m/s**2: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(DGRAV.LT.9.50D0) DGRAV=9.78032677D0*(1.D0+0.001931851353D0* 1 DSLAT**2)/DSQRT(1.D0-DEE*DSLAT**2)-0.3086D-5*DH CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute ellipsoidal curvature radius DN in meter. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DN=DEA/DSQRT(1.D0-DEE*DSLAT**2) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute geocentric latitude DPSI in degree: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DPSI=DRO*DATAN(((DN*(1.D0-DEE)+DH)*DSLAT)/((DN+DH)*DCLAT)) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute geocentric radius DR1 in meter: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DR1=DSQRT((DN+DH)**2*DCLAT**2+(DN*(1.D0-DEE)+DH)**2*DSLAT**2) IF(IPRINT.GT.0) WRITE(IUN6,7001) DLAT,DPSI,DLON,DH,DGRAV,DR1,IC, 1 DAZ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Definition of astronomical parameters I.A.U. 1984. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DPAR=3422.448D0 DMAS=1.D0/0.01230002D0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute DOODSON's constant DDC in m**2/s**2: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DR0=DEA*(1.D0-DEE/6.D0-5.D0*DEE**2/72.D0-DEE**3*55.D0/1296.D0) DDC=DR0**2*0.75D0*DEGM/(DEA**3*DMAS) DDC=DDC*(DPAR*DRAD/3600.D0)**3 DF=DRO*3.600D-3/DGRAV IF(IPRINT.GT.0) WRITE(IUN6,7002) DPAR,DMAS,DDC,DF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute latitude dependent elastic parameters from WAHR-DEHANT- C C ZSCHAU model: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL ETLOVE(IUN6,IPRINT,DLAT,DH) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DCPSI is cos and DSPSI is sin of geocentric latitude. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DCPSI=DCOS(DPSI*DRAD) DSPSI=DSIN(DPSI*DRAD) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute geodetic coefficients for tidal potential, stored C C provisional in array DGK. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DC2=DDC*(DR1/DR0)**2 DC3=DDC*(DR1/DR0)**3 DC4=DDC*(DR1/DR0)**4 DGK(1) =DC2*0.5D0*(1.D0-3.D0*DSPSI**2) DGK(2) =DC2*2.D0*DSPSI*DCPSI DGK(3) =DC2*DCPSI**2 DGK(4) =DC3*1.118033989D0*DSPSI*(3.D0-5.D0*DSPSI**2) DGK(5) =DC3*0.726184378D0*DCPSI*(1.D0-5.D0*DSPSI**2) DGK(6) =DC3*2.598076212D0*DSPSI*DCPSI**2 DGK(7) =DC3*DCPSI**3 DGK(8) =DC4*0.125000000D0*(3.D0-30.D0*DSPSI**2+35.D0*DSPSI**4) DGK(9) =DC4*0.473473091D0*2.D0*DSPSI*DCPSI*(3.D0-7.D0*DSPSI**2) DGK(10)=DC4*0.777777778D0*DCPSI**2*(1.D0-7.D0*DSPSI**2) DGK(11)=DC4*3.079201436D0*DSPSI*DCPSI**3 DGK(12)=DC4*DCPSI**4 DO 10 I=1,12 10 DPK(I)=0.D0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute geodetic coefficients for tidal acceleration vector C C orientated to sperical coordinate system, stored in DGX (north), C C DGY (east), DGZ (radially upwards), dimensions are nm/s**2. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DC2=DDC/DR0*(DR1/DR0)*1.D9 DC3=DDC/DR0*(DR1/DR0)**2*1.D9 DC4=DDC/DR0*(DR1/DR0)**3*1.D9 DGX(1) =-DC2*3.D0*DSPSI*DCPSI DGX(2) = DC2*2.D0*(DCPSI**2-DSPSI**2) DGX(3) =-DC2*2.D0*DSPSI*DCPSI DGX(4) = DC3*1.118033989D0*DCPSI*(3.D0-15.D0*DSPSI**2) DGX(5) = DC3*0.726184378D0*DSPSI*(4.D0-15.D0*DCPSI**2) DGX(6) = DC3*2.598076212D0*DCPSI*(1.D0-3.D0*DSPSI**2) DGX(7) =-DC3*3.D0*DSPSI*DCPSI**2 DGX(8) =-DC4*0.125000000D0*DSPSI*DCPSI*(60.D0-140.D0*DSPSI**2) DGX(9) = DC4*0.473473091D0*(6.D0-54.D0*DSPSI**2+56.D0*DSPSI**4) DGX(10)=-DC4*0.777777778D0*DCPSI*(16.D0*DSPSI-28.D0*DSPSI**3) DGX(11)= DC4*3.079201436D0*DCPSI**2*(4.D0*DCPSI**2-3.D0) DGX(12)=-DC4*4.D0*DCPSI**3*DSPSI C DGY(1) = 0.D0 DGY(2) = 2.D0*DC2*DSPSI DGY(3) = 2.D0*DC2*DCPSI DGY(4) = 0.D0 DGY(5) = DC3*0.726184378D0*(1.D0-5.D0*DSPSI**2) DGY(6) = 2.D0*DC3*2.598076212D0*DSPSI*DCPSI DGY(7) = 3.D0*DC3*DCPSI**2 DGY(8) = 0.D0 DGY(9) = 2.D0*DC4*0.473473091D0*DSPSI*(3.D0-7.D0*DSPSI**2) DGY(10)= 2.D0*DC4*0.777777778D0*DCPSI*(1.D0-7.D0*DSPSI**2) DGY(11)= 3.D0*DC4*3.079201436D0*DSPSI*DCPSI**2 DGY(12)= 4.D0*DC4*DCPSI**3 C DGZ(1) = DC2*(1.D0-3.D0*DSPSI**2) DGZ(2) = 4.D0*DC2*DSPSI*DCPSI DGZ(3) = 2.D0*DC2*DCPSI**2 DGZ(4) = 3.D0*DC3*1.118033989D0*DSPSI*(3.D0-5.D0*DSPSI**2) DGZ(5) = 3.D0*DC3*0.726184378D0*DCPSI*(1.D0-5.D0*DSPSI**2) DGZ(6) = 3.D0*DC3*2.598076212D0*DSPSI*DCPSI**2 DGZ(7) = 3.D0*DC3*DCPSI**3 DGZ(8) = 4.D0*DC4*0.125000000D0*(3.D0-30.D0*DSPSI**2+35.D0* 1 DSPSI**4) DGZ(9) = 8.D0*DC4*0.473473091D0*DSPSI*DCPSI*(3.D0-7.D0*DSPSI**2) DGZ(10)= 4.D0*DC4*0.777777778D0*DCPSI**2*(1.D0-7.D0*DSPSI**2) DGZ(11)= 4.D0*DC4*3.079201436D0*DSPSI*DCPSI**3 DGZ(12)= 4.D0*DC4*DCPSI**4 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute geodetic coefficients for tidal acceleration vector C C orientated to ellipsoidal coordinate system stored in C C DGX (north), DGY (east) and DGZ (upwards), all in nm/s**2. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DCDLAT=DCLAT*DCPSI+DSLAT*DSPSI DSDLAT=DSLAT*DCPSI-DCLAT*DSPSI DO 50 I=1,12 DUMMY =DCDLAT*DGX(I)-DSDLAT*DGZ(I) DGZ(I)=DSDLAT*DGX(I)+DCDLAT*DGZ(I) DGX(I)=DUMMY 50 CONTINUE IC2=IC+2 DCAZ=DCOS(DAZ*DRAD) DSAZ=DSIN(DAZ*DRAD) GOTO(100,200,300,400,500,600,700,800,900,1000,1100),IC2 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IC=-1, compute geodetic coefficients for tidal potential. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 100 CONTINUE GOTO 2000 200 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IC=0, compute geodetic coefficients for vertical component C C (gravity tide). C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 210 I=1,12 DGK(I)=DGZ(I) 210 DPK(I)=180.0D0 GOTO 2000 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IC=1, compute geodetic coefficients for horizontal component C C (tidal tilt) in azimuth DAZ, in mas. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 300 CONTINUE DO 310 I=1,12 DGK(I)=DSQRT((DGX(I)*DCAZ)**2+(DGY(I)*DSAZ)**2)*DF DPK(I)=0.D0 IF(DGX(I)*DCAZ.EQ.0.D0.AND.DGY(I)*DSAZ.EQ.0.D0) GOTO 310 DPK(I)=DRO*DATAN2(DGY(I)*DSAZ,DGX(I)*DCAZ) 310 CONTINUE GOTO 2000 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IC=2, compute geodetic coefficients for vertical displacement C C in mm. C C Attention: this component has never been tested !!! C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 400 CONTINUE DFAK=1.D3/DGRAV DO 410 I=1,12 DGK(I)=DGK(I)*DHLAT(I)*DFAK 410 DPK(I)=0.0D0 WRITE(IUN6,*) '*****The component',IC,' has never been tested !!!' WRITE(*,*) '*****The component',IC,' has never been tested !!!' GOTO 2000 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IC=3, compute geodetic coefficients for horizontal displacement C C in azimuth DAZ in mm. C C Attention: this component has never been tested !!! C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 500 CONTINUE DFAK=1.D3*DR1/DGRAV DO 510 I=1,12 DGK(I)=DSQRT((DGX(I)*DCAZ)**2+(DGY(I)*DSAZ)**2)*DLLAT(I)*DFAK DPK(I)=0.D0 IF(DGX(I)*DCAZ.EQ.0.D0.AND.DGY(I)*DSAZ.EQ.0.D0) GOTO 510 DPK(I)=DRO*DATAN2(DGY(I)*DSAZ,DGX(I)*DCAZ) 510 CONTINUE WRITE(IUN6,*) '*****The component',IC,' has never been tested !!!' WRITE(*,*) '*****The component',IC,' has never been tested !!!' GOTO 2000 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IC=4, compute geodetic coefficients for vertical strain at the C C Earth's deformed surface in 10**-9 units = nstr. C C We use a spherical approximation for the vertical strain, C C i.e. eps(rr) , and a POISSON ratio of 0.25 (see ZUERN and C C WILHELM 1984, p. 282). C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 600 CONTINUE DPOISS=0.25D0 DFAK=1.D9*DPOISS/(DPOISS-1.D0) DO 610 I=1,3 610 DGK(I)=DGK(I)*DFAK*(2.D0*DHLAT(I)-2.D0*3.D0*DLLAT(I))/(DGRAV*DR1) DO 620 I=4,7 620 DGK(I)=DGK(I)*DFAK*(2.D0*DHLAT(I)-3.D0*4.D0*DLLAT(I))/(DGRAV*DR1) DO 630 I=8,12 630 DGK(I)=DGK(I)*DFAK*(2.D0*DHLAT(I)-4.D0*5.D0*DLLAT(I))/(DGRAV*DR1) DO 640 I=1,12 640 DPK(I)=0.0D0 GOTO 2000 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IC=5, compute geodetic coefficients for horizontal strain C C in azimuth DAZ, in 10**-9 units. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 700 CONTINUE DTHETA=(90.D0-DPSI)*DRAD DAZR=(DAZ+180.D0)*DRAD DCAZ =DCOS(DAZR) DSAZ =DSIN(DAZR) DSAZ2=DSIN(2.D0*DAZR) DCSTS=-0.5D0*DSIN(2.D0*DAZR) DCT=DSPSI DST=DCPSI DCT2=DCT*DCT DST2=DST*DST DCC2=DCOS(2.D0*DPSI*DRAD) DC2T=-DCC2 DCOTT =1.D0/DTAN(DTHETA) DCOTT2=1.D0/DTAN(2.D0*DTHETA) DFAK=1.D9/(DR1*DGRAV) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Real part is stored in DGX, imaginary part is stored in DGY. C C Formulas were given by Dr. W. Zuern, BFO Schiltach (personal C C communication) and tested against horizontal strain computed C C (with lower precision) by program ETIDEL. C C Results agreed to 0.3 % and 0.1 degree for most of the waves, C C except for 2N2 and L2 (deviation of 3 %). C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DGX(1)=(DHLAT(1)-(6.D0*DLLAT(1)*DC2T)/(3.D0*DCT2-1.D0))*DCAZ**2 1 +(DHLAT(1)-(6.D0*DLLAT(1)*DCT2)/(3.D0*DCT2-1.D0))*DSAZ**2 DGY(1)=0.D0 DGX(2)=(DHLAT(2)-4.D0*DLLAT(2))*DCAZ**2+(DHLAT(2)-DLLAT(2)/DST2 1 +2.D0*DLLAT(2)*DCOTT*DCOTT2)*DSAZ**2 DGY(2)=2.D0*DLLAT(2)*(2.D0*DCOTT2-DCOTT)*DCSTS/DST DGX(3)=(DHLAT(3)+2.D0*DLLAT(3)*(DCOTT*DCOTT-1.D0))*DCAZ**2 1 +(DHLAT(3)-4.D0*DLLAT(3)/DST2+2.D0*DLLAT(3)*DCOTT*DCOTT)*DSAZ**2 DGY(3)=4.D0*DLLAT(3)*DCOTT*DCSTS/DST DGX(4)=(DHLAT(4)+DLLAT(4)*(33.D0-45.D0*DCT2)/(5.D0*DCT2-3.D0))* 1 DCAZ**2+(DHLAT(4)-DLLAT(4)*(1.D0+10.D0*DCT2/(5.D0*DCT2-3.D0)))* 2 DSAZ**2 DGY(4)=0.D0 DGX(5)=(DHLAT(5)-DLLAT(5)*(1.D0+10.D0*(1.D0-4.D0*DCT2)/ 1 (1.D0-5.D0*DCT2)))*DCAZ**2+(DHLAT(5)+DLLAT(5)* 2 (DCOTT*DCOTT-1.D0/DST2-10.D0*DCT2/(5.D0*DCT2-1.D0)))*DSAZ**2 DGY(5)=-20.D0*DLLAT(5)*DCT*DCSTS/(5.D0*DCT2-1.D0) DGX(6)=(DHLAT(6)+DLLAT(6)*(2.D0*DCOTT*DCOTT-7.D0))*DCAZ**2 1 +(DHLAT(6)+DLLAT(6)*(2.D0*DCOTT*DCOTT-1.D0-4.D0/DST2))*DSAZ**2 DGY(6)=-4.D0*DLLAT(6)*(DCOTT-1.D0/DCOTT)*DCSTS/DST DGX(7)=(DHLAT(7)+DLLAT(7)*(6.D0*DCOTT*DCOTT-3.D0))*DCAZ**2 1 +(DHLAT(7)+DLLAT(7)*(3.D0*DCOTT*DCOTT-9.D0/DST2))*DSAZ**2 DGY(7)=12.D0*DLLAT(7)*DCOTT*DCSTS/DST DGX(8)=(DHLAT(8)-4.D0*DLLAT(8)*(4.D0-3.D0*(5.D0*DCT2-1.D0)/ 1 (35.D0*DCT2*DCT2-30.D0*DCT2+3.D0)))*DCAZ**2+ 2 (DHLAT(8)-4.D0*DLLAT(8)*(1.D0+3.D0*(5.D0*DCT2-1.D0)/ 3 (35.D0*DCT2*DCT2-30.D0*DCT2+3.D0)))*DSAZ**2 DGY(8)=0.D0 DGX(9)= (DHLAT(9)-2.D0*DLLAT(9)*(8.D0-3.D0/(7.D0*DCT2-3.D0)))* 1 DCAZ**2+(DHLAT(9)-2.D0*DLLAT(9)*(2.D0+3.D0/(7.D0*DCT2-3.D0)))* 2 DSAZ**2 DGY(9)=DLLAT(9)*3.D0/DCT*(1.D0+2.D0/(7.D0*DCT2-3.D0))*DSAZ2 DGX(10)=(DHLAT(10)-4.D0*DLLAT(10)*(4.D0+3.D0*DCT2/ 1 (7.D0*DCT2**2-8.D0*DCT2+1.D0)))*DCAZ**2 2 +(DHLAT(10)-4.D0*DLLAT(10)*(1.D0-3.D0*DCT2/ 2 (7.D0*DCT2**2-8.D0*DCT2+1.D0)))*DSAZ**2 DGY(10)=-DLLAT(10)*6.D0*DCT/DST**2*(1.D0-4.D0/(7.D0*DCT2-1.D0))* 1 DSAZ2 DGX(11)=(DHLAT(11)-2.D0*DLLAT(11)*(8.D0-3.D0/DST2))*DCAZ**2 1 +(DHLAT(11)-2.D0*DLLAT(11)*(2.D0+3.D0/DST2))*DSAZ**2 DGY(11)= DLLAT(11)*3.D0/DCT*(3.D0-2.D0/DST2)*DSAZ2 DGX(12)=(DHLAT(12)-4.D0*DLLAT(12)*(4.D0-3.D0/DST2))*DCAZ**2 1 +(DHLAT(12)-4.D0*DLLAT(12)*(1.D0+3.D0/DST2))*DSAZ**2 DGY(12)= DLLAT(12)*12.D0*DCT/DST2*DSAZ2 DO 710 I=1,12 DGK(I)=DGK(I)*DSQRT(DGX(I)**2+DGY(I)**2)*DFAK 710 DPK(I)=DPK(I)+DATAN2(DGY(I),DGX(I))*DRO GOTO 2000 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IC=6, compute geodetic coefficients for areal strain C C in 10**-9 units = nstr. C C We use a spherical approximation for the aereal strain, C C i.e. eps(t,t) + eps(l,l), (see ZUERN and WILHELM 1984, C C p. 282). C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 800 CONTINUE DO 810 I=1,3 810 DGK(I)=DGK(I)*(2.D0*DHLAT(I)-2.D0*3.D0*DLLAT(I))/(DGRAV*DR1)*1.D9 DO 820 I=4,7 820 DGK(I)=DGK(I)*(2.D0*DHLAT(I)-3.D0*4.D0*DLLAT(I))/(DGRAV*DR1)*1.D9 DO 830 I=8,12 830 DGK(I)=DGK(I)*(2.D0*DHLAT(I)-4.D0*5.D0*DLLAT(I))/(DGRAV*DR1)*1.D9 DO 840 I=1,12 840 DPK(I)=0.0D0 GOTO 2000 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IC=7, compute geodetic coefficients for shear tidal strain C C at the Earth's deformed surface in 10**-9 units = nstr. C C We use a spherical approximation, i.e. eps(t,l) C C Attention: this component has never been tested !!!! C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 900 CONTINUE DTHETA=(90.D0-DPSI)*DRAD DAZR=(DAZ+180.D0)*DRAD DCAZ =DCOS(DAZR) DSAZ =DSIN(DAZR) DSAZ2=DSIN(2.D0*DAZR) DCSTS=-0.5D0*DSIN(2.D0*DAZR) DCT=DSPSI DST=DCPSI DCT2=DCT*DCT DST2=DST*DST DCC2=DCOS(2.D0*DPSI*DRAD) DC2T=-DCC2 DCOTT =1.D0/DTAN(DTHETA) DCOTT2=1.D0/DTAN(2.D0*DTHETA) DFAK=1.D9/(DR1*DGRAV) DGY(1)=0.D0 DGY(2)=2.D0*DLLAT(2)*(2.D0*DCOTT2-DCOTT)*DCSTS/DST DGY(3)=4.D0*DLLAT(3)*DCOTT*DCSTS/DST DGY(4)=0.D0 DGY(5)=-20.D0*DLLAT(5)*DCT*DCSTS/(5.D0*DCT2-1.D0) DGY(6)=-4.D0*DLLAT(6)*(DCOTT-1.D0/DCOTT)*DCSTS/DST DGY(7)=12.D0*DLLAT(7)*DCOTT*DCSTS/DST DGY(8)=0.D0 DGY(9)=DLLAT(9)*3.D0/DCT*(1.D0+2.D0/(7.D0*DCT2-3.D0))*DSAZ2 DGY(10)=-DLLAT(10)*6.D0*DCT/DST**2*(1.D0-4.D0/(7.D0*DCT2-1.D0))* 1 DSAZ2 DGY(11)=DLLAT(11)*3.D0/DCT*(3.D0-2.D0/DST2)*DSAZ2 DGY(12)=DLLAT(12)*12.D0*DCT/DST2*DSAZ2 DO 910 I=1,12 DGK(I)=DGK(I)*DGY(I)*DFAK 910 DPK(I)=0.D0 WRITE(IUN6,*) ' ***** The shear strain has never been tested !!!' WRITE(*,*) ' ***** The shear strain has never been tested !!!' GOTO 2000 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IC=8, compute geodetic coefficients for volume strain C C at the Earth's deformed surface in 10**-9 units = nstr. C C We use a spherical approximation, i.e. eps(t,t)+eps(l,l). C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1000 CONTINUE DPOISS=0.25D0 DFAK=1.D9*(1.D0-2.D0*DPOISS)/(1.D0-DPOISS) DO 1010 I=1,3 1010 DGK(I)=DGK(I)*DFAK*(2.D0*DHLAT(I)-2.D0*3.D0*DLLAT(I))/(DGRAV*DR1) DO 1020 I=4,7 1020 DGK(I)=DGK(I)*DFAK*(2.D0*DHLAT(I)-3.D0*4.D0*DLLAT(I))/(DGRAV*DR1) DO 1030 I=8,12 1030 DGK(I)=DGK(I)*DFAK*(2.D0*DHLAT(I)-4.D0*5.D0*DLLAT(I))/(DGRAV*DR1) DO 1040 I=1,12 1040 DPK(I)=0.0D0 GOTO 2000 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IC=9, compute geodetic coefficients for ocean tides in mm: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1100 CONTINUE DFAK=1.D3/DGRAV DO 1110 I=1,12 DGK(I)=DGK(I)*DFAK 1110 DPK(I)=0.0D0 2000 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Common loop for all components. C C For negative geodetic coefficients, use absolute value C C and add 180 deg to the phase. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 2010 I=1,12 IF(DGK(I).GE.0.D0) GOTO 2010 DGK(I)=-DGK(I) DPK(I)=DPK(I)+180.D0 IF(DPK(I).LT.0.D0) DPK(I)=DPK(I)+360.D0 2010 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Print geodetic coefficients: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(IPRINT.EQ.0) RETURN WRITE(IUN6,7003) IC,DAZ,(DGK(I),CUNIT(IC2),DPK(I),I=1,12) 5000 CONTINUE IF(IPRINT.GT.0) WRITE(IUN6,7005) RETURN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Format statements: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 7000 FORMAT(' Routine ETGCOF, version 930706 FTN 77.'// 1' Computation of geodetic coefficients'// 3' Parameters of Geodetic Reference System 1980:'/ 4' Major semi axis ',F12.0,' m'/ 5' 1. excentricity ',F12.8/ 6' Geocentric gravitational constant',F12.1,' 10**9 m**3/s**2'/) 7001 FORMAT(' Station paremeters:'// 1' Latitude ',F12.6,' deg'/ 2' Geocentric latitude ',F12.6,' deg'/ 3' Longitude ',F12.6,' deg'/ 4' Height ',F12.3,' m'/ 5' Gravity ',F12.6,' m/s**2'/ 6' Geocentric radius ',F12.3,' m'/ 7' Component of observations ',I12/ 8' Azimuth from north direction ',F12.6,' deg'//) 7002 FORMAT(' Astronomic constants International Astronomic Union', 1' 1984 m.kg.s-system'/ 2' Moons mean sine parallax ',F12.6/ 3' Mass relation Earth/Moon ',F12.6// 4' DOODSONs constant ',F12.6,' m**2/s**2'/ 6' F ',F12.6,' mas/(nm/s**2)'//) 7003 FORMAT(//' Geodetic coefficients and phases for component',I4/ 1' azimuth:',F12.6,' degree'// 2' GC 2,0',F14.5,2X,A8,2X,F14.6,' deg'/ 3' GC 2,1',F14.5,2X,A8,2X,F14.6,' deg'/ 4' GC 2,2',F14.5,2X,A8,2X,F14.6,' deg'/ 5' GC 3,0',F14.5,2X,A8,2X,F14.6,' deg'/ 6' GC 3,1',F14.5,2X,A8,2X,F14.6,' deg'/ 7' GC 3,2',F14.5,2X,A8,2X,F14.6,' deg'/ 8' GC 3,3',F14.5,2X,A8,2X,F14.6,' deg'/ 9' GC 4,0',F14.5,2X,A8,2X,F14.6,' deg'/ *' GC 4,1',F14.5,2X,A8,2X,F14.6,' deg'/ 1' GC 4,2',F14.5,2X,A8,2X,F14.6,' deg'/ 2' GC 4,3',F14.5,2X,A8,2X,F14.6,' deg'/ 3' GC 4,4',F14.5,2X,A8,2X,F14.6,' deg'//) 7005 FORMAT(///' ***** Routine ETGCOF finished the execution.'/) END C SUBROUTINE ETGREI(ITY,ITM,ITD,ITH) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Routine ETGREI, version 901230 FORTRAN 77. C C C C The routine ETGREI computes the GREGORIAN date from year, C C month, day and hour, where the hour may exceed 24. C C C C Input/output parameter description: C C ----------------------------------- C C C C All following parameters are input and output parameters, C C which measn that the parameters may be changed during C C the execution of routine ETGREI. C C C C ITY... year in INTEGER form, E.G. 1971 C C ITM... month in INTEGER form, E.G. 1 = January. C C ITD... day in INTEGER form. C C ITH... hour in INTEGER form (UTC). ITH may exceed 24. C C C C Routine creation: 710523 by H.-G. Wenzel. C C Geodaetisches Institut, C C Universitaet Karlsruhe, C C Englerstr. 7, C C D-76128 KARLSRUHE 1, C C Germany. C C Tel: 0049-721-6082307, C C FAX: 0049-721-694552. C C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de C C Last modification: 901230 by H.-G. Wenzel. C C**********************************************************************C IMPLICIT INTEGER*4 (I-N) INTEGER*2 ID1(12),ID2(12) SAVE ID1,ID2 DATA ID1/31,28,31,30,31,30,31,31,30,31,30,31/ DATA ID2/31,29,31,30,31,30,31,31,30,31,30,31/ IH=ITH/24 ITH=ITH-24*IH ITD=ITD+IH 50 L=ITY/4 L=4*L IF(L.EQ.ITY) GOTO 300 100 IF(ITD.LE.ID1(ITM)) GOTO 555 ITD=ITD-ID1(ITM) ITM=ITM+1 IF(ITM.EQ.13) GOTO 200 GOTO 100 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Next year: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 200 ITM=1 ITY=ITY+1 GOTO 50 300 IF(ITD.LE.ID2(ITM)) GOTO 555 ITD=ITD-ID2(ITM) ITM=ITM+1 IF(ITM.EQ.13) GOTO 200 GOTO 300 555 RETURN END C SUBROUTINE ETINPD(IUN5,IUN6,IUN20,IPRINT,NC,NB,NREC,NERR) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Routine ETINPD, version 940117 FORTRAN 77. C C C C The routine ETINPD reads observations from input iunit IUN5 and C C stores them on direct access unit IUN20. C C C C Input parameter description: C C ---------------------------- C C C C IUN5... formatted input data unit. C C IUN6... formatted printer unit. C C IUN20... unformatted direct access unit. C C IPRINT... printout parameter. For IPRINT = 0, noting will be C C written to printer unit IUN6. C C NC... number of data channels (including Earth tide data C C and meteorological data). NC is restricted to be C C or equal to 6. C C C C Output parameter description: C C ----------------------------- C C C C NB... number of blocks (uninterrupted parts of data). C C NB is restricted to be less or equal to 300. C C NREC... total number of data records stored on direct C C access unit IUN20. NREC is not restricted. C C NERR... number of sequence errors ocuured during input of C C data. C C C C Description of COMMON (BLOCKR/: C C ------------------------------- C C C C IDATA... array(1:300) of start dates of the blocks. C C IDATE... array(1:300) of end dates of the blocks. C C C C Program creation: 930629 by H.-G. Wenzel, C C Geodaetisches Institut, C C Universitaet Karlsruhe, C C Englerstr. 7, C C D-76128 KARLSRUHE 1, C C Germany. C C Tel.: 0721-6082301. C C FAX: 0721-694552. C C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de C C Last Modification: 940117 by H.-G.Wenzel. C C**********************************************************************C IMPLICIT DOUBLE PRECISION (D) CHARACTER C8888*8 DIMENSION DCIN(6),DOFFS(6) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C The following dimension statements are concerning the number of C C blocks of data without interruption, which is restricted to 300 C C in the current program version. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC INTEGER IRECA(300),IRECE(300),IDATA(300),IDATE(300),ITIMA(300), 1 ITIME(300),IOB(300),NBIAS(300) DOUBLE PRECISION DSAPR(300),DSAPO(300),DTLAG(300) CHARACTER CINSTR(300)*10 COMMON /BLOCKR/ IRECA,IRECE,IDATA,ITIMA,IDATE,ITIME,IOB,NBIAS, 1 DSAPR,DSAPO,DTLAG COMMON /BLOCKC/ CINSTR DATA C8888/'88888888'/ NB=0 IREC=0 NERR=0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C New block: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 100 CONTINUE NB=NB+1 IREC=IREC+1 READ(IUN5,7002) CINSTR(NB),DCAL,DSAPR(NB),DTLAG(NB),NBIAS(NB) IF(CINSTR(NB).EQ.C8888) GOTO 1000 IF(KFILT.GT.0) NBIAS(NB)=0 IRECA(NB)=IREC IRECE(NB)=IREC READ(IUN5,7000) (DOFFS(J),J=1,NC) READ(IUN5,7001) IDAT,ITIM,(DCIN(J),J=1,NC) DO 110 J=1,NC 110 DCIN(J)=DCIN(J)+DOFFS(J) IF(IPRINT.EQ.0) GOTO 120 WRITE(*,7003) IDAT,ITIM,(DCIN(J),J=1,NC) WRITE(IUN6,7003) IDAT,ITIM,(DCIN(J),J=1,NC) 120 DCIN(1)=DCIN(1)*DCAL WRITE(IUN20,REC=IREC) IDAT,ITIM,(DCIN(J),J=1,NC) IOB(NB)=1 IDATA(NB)=IDAT ITIMA(NB)=ITIM IDATE(NB)=IDAT ITIME(NB)=ITIM ITIMN=ITIM+10000 IDATN=IDAT IF(ITIMN.LT.230000) GOTO 130 CALL DATUM(IDATN,ITIMN,ITYN,ITMN,ITDN,ITHN) CALL ETGREI(ITYN,ITMN,ITDN,ITHN) IDATN=ITDN+100*ITMN+10000*ITYN ITIMN=10000*ITHN 130 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Read rest of the block: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 200 CONTINUE READ(IUN5,7001) IDAT,ITIM,(DCIN(J),J=1,NC) IF(IDAT.EQ.99999999) GOTO 100 IF(IDAT.NE.77777777) GOTO 220 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Step, update offsets: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 210 JC=1,NC 210 DOFFS(JC)=DOFFS(JC)+DCIN(JC) GOTO 200 220 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Check the sequence of data: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(ITIM.EQ.ITIMN.AND.IDAT.EQ.IDATN) GOTO 230 WRITE(*,7004) IDAT,ITIM WRITE(IUN6,7004) IDAT,ITIM NERR=NERR+1 230 CONTINUE ITIMN=ITIMN+10000 IF(ITIMN.LT.230000) GOTO 240 CALL DATUM(IDATN,ITIMN,ITYN,ITMN,ITDN,ITHN) CALL ETGREI(ITYN,ITMN,ITDN,ITHN) IDATN=ITDN+100*ITMN+10000*ITYN ITIMN=10000*ITHN 240 CONTINUE DO 250 J=1,NC 250 DCIN(J)=DCIN(J)+DOFFS(J) IF(IPRINT.EQ.0) GOTO 260 WRITE(*,7003) IDAT,ITIM,(DCIN(J),J=1,NC) WRITE(IUN6,7003) IDAT,ITIM,(DCIN(J),J=1,NC) 260 DCIN(1)=DCIN(1)*DCAL IREC=IREC+1 WRITE(IUN20,REC=IREC) IDAT,ITIM,(DCIN(J),J=1,NC) IOB(NB)=IOB(NB)+1 IRECE(NB)=IREC IDATE(NB)=IDAT ITIME(NB)=ITIM GOTO 200 1000 CONTINUE NB=NB-1 NREC=IREC-1 RETURN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Format statements: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 7000 FORMAT(15X,6F10.3) 7001 FORMAT(I8,1X,I6,6F10.3) 7002 FORMAT(A10,5X,2F10.4,F10.3,I10) 7003 FORMAT(1X,I8,1X,I6,6F10.3) 7004 FORMAT(' ***** Error of sequence at ',I9,1X,I6) END C SUBROUTINE ETJULD(IUN6,ITY,ITM,ITD,DTH,DJULD) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Routine ETJULD, version 930503 FORTRAN 77. C C C C The routine ETJULD computes the Julian date from year, C C month, days and hours. C C C C Input parameter description: C C ---------------------------- C C C C IUN6... formatted printer unit. C C ITY... year in integer form, e.g. 1971 C C ITM... month in INTEGER form, e.g. 1 = JANUARY. C C ITD... day in integer form. C C DTH... double precision hour in UTC. C C C C Output parameter description: C C ----------------------------- C C C C DJULD... double precision Julian date, e.g. for C C 01. January 1988, 0.00 H UTC is DJULD=2447161.500 C C 01. February 1988, 0.00 H UTC is DJULD=2447192.500 C C 29. February 1988, 0.00 H UTC is DJULD=2447220.500 C C 01. March 1988, 0.00 H UTC is DJULD=2447221.500 C C C C Execution time: C C --------------- C C C C About 0.0062 sec CPU time for 1000 calls of routine ETJULD on C C 80486 DX2 66 Mhz processor under operation system MS-DOS using C C the LAHEY FORTRAN-compiler. C C C C Numerical accuracy: C C ------------------- C C C C The result under operation system MS-DOS is accurate to C C 0.000 035 s which is equivalent to 4*E-11 days. C C C C Program creation: 710523 by H.-G. Wenzel, C C Geodaetisches Institut, C C Universitaet Karlsruhe, C C Englerstr. 7, C C D-76128 KARLSRUHE 1, C C Germany. C C Tel.: 0721-6082301. C C FAX: 0721-694552. C C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de C C Last Modification: 930503 by H.-G.Wenzel. C C**********************************************************************C IMPLICIT DOUBLE PRECISION (D) IMPLICIT INTEGER*4 (I-N) INTEGER*2 ITDY(12) SAVE ITDY DATA ITDY/0,31,59,90,120,151,181,212,243,273,304,334/ DJULD=2415019.5D0 ITYR=ITY-1900 ILEAP=ITYR/4 ITA=ITYR*365+ILEAP IF(ITM.LT.1) GOTO 5000 IF(ITM.GT.12) GOTO 5010 DJULD=DJULD+DBLE(ITA+ITDY(ITM)+ITD)+DTH/24.D0 ITYL=ILEAP*4 IF(ITYR.NE.ITYL) RETURN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Leap year: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(ITM.GT.2) RETURN DJULD=DJULD-1.D0 RETURN 5000 WRITE(IUN6,7050) ITY,ITM,ITD,DTH WRITE(*,7050) ITY,ITM,ITD,DTH STOP 5010 WRITE(IUN6,7051) ITY,ITM,ITD,DTH WRITE(*,7051) ITY,ITM,ITD,DTH STOP CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Format statements: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 7050 FORMAT(/' *****Error in routine ETJULD, version 930503 FTN77.'/ 1' *****Month is less 1:',2X,3I4,F12.3/ 2' *****Routine ETJULD stops the execution.'/) 7051 FORMAT(/' *****Error in routine ETJULD, version 930503 FTN77.'/ 1' *****Month is greater 12:',2X,3I4,F12.3/ 2' *****Routine ETJULD stops the execution.'/) END C SUBROUTINE ETLOVE(IUN6,IPRINT,DLAT,DELV) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Routine ETLOVE, version 930706 FORTRAN 77. C C C C The routine computes latitude dependent LOVE-numbers DH, DK, C C SHIDA-numbers DL, gravimeter factors DG and tilt factors DT C C using the WAHR-DEHANT-ZSCHAU model. C C C C Body tide amplitude factors for WAHR-DEHANT-ZSCHAU model. C C The NDFW resonance is approximated by C C C C G(RES) = GLAT - GR*(DOM - DOM0)/(DOMR - DOM). C C C C similar equations hold for the other parameters. C C C C Gravimetric amplitude factors, LOVE numbers h and k for degree C C 0...3 have been taken from DEHANT 1987, Table 7, 8 and 9 C C for an elliptical, uniformly rotating, oceanless Earth with C C liquid outer core and inelastic mantle (PREM Earth model with C C inelastic mantle from ZSCHAU) and for the fourth degree from C C DEHANT et. al 1989, Table 6). The resonance factors GR have C C been computed to fit the difference between body tide amplitude C C factors at O1 and PSI1 from DEHANT 1987, PREM model with C C elastic mantle (Table 1...3). The NDFW resonance frequency is C C 15.073729 degree per hour = 1.004915267 CPD UT, taken from C C WAHR 1981 (because it is not given in DEHANT's papers). C C C C Input parameter description: C C ---------------------------- C C C C IUN6... formatted line printer unit. C C IPRINT... printout parameter. For IPRINT=1, the computed C C LOVE- and SHIDA- number s will be printed. C C DLAT... ellipsoidal latitude in degree. C C DELV... ellipsoidal height in meter. C C C C Description of COMMON /LOVE/: C C ----------------------------- C C C C DOM0... frequency of O1 in degree per hour. C C DOMR... frequency of the FCN eigenfrequency in degree per C C hour. C C DGLAT... array(1..12) containing the gravimetric factors at C C latitude DLAT. C C DGR... resonance factor for gravimetric factors. C C DHLAT array(1..12) containing the LOVE-numbers h at C C latitude DLAT. C C DHR... resonance factor for the LOVE-number h(2,1). C C DKLAT array(1..12) containing the LOVE-numbers k at C C latitude DLAT. C C DKR... resonance factor for the LOVE-number k(2,1). C C DLLAT... array(1..12) containing the SHIDA-numbers l at C C latitude DLAT. C C DLR... resonance factor for the SHIDA-number l(2,1). C C DTLAT... array(1..12) containing the tilt factors at C C latitude DLAT. C C C C Reference: C C ---------- C C C C C C DEHANT, V. 1987: Tidal Parameters for an Inelastic Earth. C C Physics of the Earth and Planetary Interiors, 49, 97-116, C C 1987. C C C C WAHR, J.M. 1981: Body tides on an elliptical, rotating, elastic C C and oceanless earth. Geophysical Journal of the Royal C C astronomical Society, vol. 64, 677-703, 1981. C C C C ZSCHAU, J. and R. WANG 1987: Imperfect elasticity in the Earth's C C mantle. Implications for Earth tides and long period C C deformations. Proceedings of the 9th International Sym- C C posium on Earth Tides, New York 1987, pp. 605-629, editor C C J.T. KUO, Schweizerbartsche Verlagsbuchhandlung, Stuttgart C C 1987. C C C C C C Routine creation: 930703 by H.-G. Wenzel. C C Geodaetisches Institut, C C Universitaet Karlsruhe, C C Englerstr. 7, C C D-76128 KARLSRUHE 1, C C Germany. C C Tel: 0049-721-6082307, C C FAX: 0049-721-694552. C C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de C C Last modification: 930706 by H.-G. Wenzel. C C**********************************************************************C IMPLICIT DOUBLE PRECISION (D) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C The following DIMENSION statement is concerning the elastic C C Earth model for the different degree and order constituents. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DOUBLE PRECISION DG0(12),DGP(12),DGM(12) DOUBLE PRECISION DH0(12),DHP(12),DHM(12) DOUBLE PRECISION DK0(12),DKP(12),DKM(12) DOUBLE PRECISION DL0(12),DLP(12),DLM(12) DOUBLE PRECISION DLATP(12),DLATM(12) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C COMMON /LOVE/ contains gravimeter factors, LOVE-numbers, SHIDA- C C numbers and tilt factors for degree 2...4 at latitude DLAT: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DIMENSION DGLAT(12),DHLAT(12),DKLAT(12),DLLAT(12),DTLAT(12) COMMON /LOVE/ DOM0,DOMR,DGLAT,DGR,DHLAT,DHR,DKLAT,DKR,DLLAT,DLR, 1 DTLAT,DTR CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C COMMON /CONST/: C C DPI... 3.1415.... C C DPI2... 2.D0*DPI C C DRAD... DPI/180.D0 C C DRO... 180.D0/DPI C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC COMMON /CONST/ DPI,DPI2,DRAD,DRO CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C The following DATA statements are concerning the elastic C C Earth model for the different degree and order constituents. C C The latitude dependency is not given for all constituents in C C the WAHR-DEHANT-ZSCHAU model !!!!!! C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DATA DG0/1.1576D0,1.1542D0,1.1600D0,1.0728D0,1.0728D0,1.0728D0, 1 1.0728D0,1.0363D0,1.0363D0,1.0363D0,1.0363D0,1.0363D0/ DATA DGP/-0.0016D0,-0.0018D0,-0.0010D0,0.D0,0.D0,0.D0,-0.0010D0, 1 0.D0,0.D0,0.D0,0.D0,-0.000315D0/ DATA DGM/0.0054D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 1 0.D0,0.D0/ DATA DH0/0.6165D0,0.6069D0,0.6133D0,0.2946D0,0.2946D0,0.2946D0, 1 0.2946D0,0.1807D0,0.1807D0,0.1807D0,0.1807D0,0.1807D0/ DATA DHP/0.0007D0,0.0007D0,0.0005D0,0.D0,0.D0,0.D0,0.0003D0, 1 0.D0,0.D0,0.D0,0.D0,0.00015D0/ DATA DHM/0.0018D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 1 0.D0,0.D0/ DATA DK0/0.3068D0,0.3009D0,0.3034D0,0.0942D0,0.0942D0,0.0942D0, 1 0.0942D0,0.0427D0,0.0427D0,0.0427D0,0.0427D0,0.0427D0/ DATA DKP/0.0015D0,0.0014D0,0.0009D0,0.D0,0.D0,0.D0,0.0007D0, 1 0.D0,0.D0,0.D0,0.D0,0.00066D0/ DATA DKM/-0.0004D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 1 0.D0,0.D0/ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SHIDA-numbers: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DATA DL0/ 0.0840D0,0.0841D0,0.0852D0,0.0149D0,0.0149D0,0.0149D0, 1 0.0149D0,0.0100D0,0.0100D0,0.0100D0,0.0100D0,0.0100D0/ DATA DLP/-0.002D0,-0.002D0,-0.001D0,0.0000D0,0.0000D0,0.0000D0, 1 0.0000D0,0.0000D0,0.0000D0,0.0000D0,0.0000D0,0.0000D0/ DATA DLM/ 0.0000D0,0.0000D0,0.0000D0,0.0000D0,0.0000D0,0.0000D0, 1 0.0000D0,0.0000D0,0.0000D0,0.0000D0,0.0000D0,0.0000D0/ DATA DLATP/0.D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 1 0.D0,0.D0/ DATA DLATM/0.D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0, 1 0.D0,0.D0/ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Definition of parameters of Geodetic Reference System 1980. C C DEA is major semi axis in meter. C C DEE is square of first excentricity (without dimnension). C C DEGM is geocentric gravitational constant in m*3/s**2. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DATA DEA/6378137.00D0/,DEE/6.69438002290D-3/ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Define resonance frequency and resonance factors: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DOMR=15.073729D0 DOM0=13.943036D0 DGR =-0.000625D0 DHR =-0.002505D0 DKR =-0.001261D0 DLR =0.0000781D0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DCLAT is cos and DSLAT is sin of ellipsoidal latitude. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DCLAT=DCOS(DLAT*DRAD) DSLAT=DSIN(DLAT*DRAD) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute ellipsoidal curvature radius DN in meter. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DN=DEA/DSQRT(1.D0-DEE*DSLAT**2) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute geocentric latitude DPSI in degree: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DPSI=DRO*DATAN(((DN*(1.D0-DEE)+DELV)*DSLAT)/((DN+DELV)*DCLAT)) DTHET=90.D0-DPSI DCT=DCOS(DTHET*DRAD) DCT2=DCT*DCT DLATP(1)=0.335410D0*(35.D0*DCT2*DCT2-30.D0*DCT2+3.D0)/ 1 (3.D0*DCT2-1.D0) DLATM(1) =0.894427D0/(3.D0*DCT2-1.D0) DLATP(2) =0.612372D0*(7.D0*DCT2-3.D0) DLATP(3) =0.866025D0*(7.D0*DCT2-1.D0) DLATP(7) =0.829156D0*(9.D0*DCT2-1.D0) DLATP(12)=0.806226D0*(11.D0*DCT2-1.D0) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute latitude dependent gravimeter factors DG: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 110 I=1,12 110 DGLAT(I)=DG0(I)+DGP(I)*DLATP(I)+DGM(I)*DLATM(I) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute latitude dependent LOVE-numbers DH (for vertical C C displacement): C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 120 I=1,12 120 DHLAT(I)=DH0(I)+DHP(I)*DLATP(I)+DHM(I)*DLATM(I) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute latitude dependent LOVE-numbers DK: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 130 I=1,12 130 DKLAT(I)=DK0(I)+DKP(I)*DLATP(I)+DKM(I)*DLATM(I) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute latitude dependent SHIDA-numbers DL: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 140 I=1,12 140 DLLAT(I)=DL0(I)+DLP(I)*DLATP(I)+DLM(I)*DLATM(I) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute latitude dependent tilt factors DT: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 150 I=1,12 DTLAT(I)=1.D0+DK0(I)-DH0(I)+DLATP(I)*(DKP(I)-DHP(I))+ 1 DLATM(I)*(DKM(I)-DHM(I)) 150 CONTINUE DTR=DKR-DHR IF(IPRINT.EQ.0) RETURN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Print out of parameters: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC WRITE(IUN6,7001) DOM0,DOMR,DGR,DHR,DKR,DLR,DTR I=0 WRITE(IUN6,7002) DLAT DO 300 L=2,4 WRITE(IUN6,7004) DO 300 M=0,L I=I+1 WRITE(IUN6,7003) L,M,DGLAT(I),DHLAT(I),DKLAT(I),DLLAT(I),DTLAT(I) 300 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Format statements: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 7001 FORMAT(/' Routine ETLOVE, version 930706 FTN77.'/ 1' Latitude dependent parameters for an elliptical, rotating,'/ 2' inelastic and oceanless Earth from WAHR-DEHANT-ZSCHAU model.'// 3' frequency of wave O1:',F10.6,' deg per hour'/ 4' resonance frequency :',F10.6,' deg per hour'// 5' resonance factor for G:',F10.6/ 6' resonance factor for h:',F10.6/ 7' resonance factor for k:',F10.6/ 8' resonance factor for l:',F10.6/ 9' resonance factor for T:',F10.6/) 7002 FORMAT(//' Latitude dependent elastic parameters'// 1' ellipsoidal latitude:',F10.4,' deg'// 2' G is gravimetric factor delta'/ 3' h is LOVE-number h'/ 4' k is LOVE-number k'/ 5' l is SHIDA-number l'/ 6' T is tilt factor gamma'// 7' degree order G h k l', 8' T') 7003 FORMAT(2I7,5F10.6) 7004 FORMAT(' ') RETURN END C SUBROUTINE ETMUTC(IUN6,IPRINT,DTUJD,DDT) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Routine ETMUTC, version 930712 FORTRAN 77. C C C C All variables with D as first character are DOUBLE PRECISION. C C C C Input parameter description: C C ---------------------------- C C C C IUN6... formatted printer unit. C C IPRINT... printout parameter. For IPRINT=0, nothing will be C C written on unit IUN6. For IPRINT=2, the tables DTX, C C DTJULD and DTY will be printed for the first call C C of routine ETMUTC. C C DTUJD... Julian date of epoch (Universal time). C C C C Output parameter description: C C ----------------------------- C C C C DDT... difference ET - UTC resp. TDT - UTC in seconds C C from 1955.5 until 1992.5 . For epochs less 1955.5, C C DDT is set to 31.59 s. For epochs exceeding 1993.5, C C DDT is set to 60.184 s. ET is Ephemeris Time. C C TDT is Terrestrial Dynamical Time. C C UTC is Universal Time Coordinated, as bradcasted by C C radio or GPS satellites. C C C C The table DTAB has to be extended, when new data are available. C C Change parameter NTAB and DIMENSIONS !!! C C C C Routine creation: W. Zuern, BFO Schiltach, Heubach 206, C C D-7620 WOLFACH (in program RIGTID). C C Last modification: 930712 by H.-G. Wenzel. C C Geodaetisches Institut, C C Universitaet Karlsruhe, C C Englerstr. 7, C C D-76128 KARLSRUHE 1, C C Germany. C C Tel.: 0721-6082307, C C FAX: 0721-694552. C C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de C C**********************************************************************C IMPLICIT DOUBLE PRECISION (D) IMPLICIT INTEGER*4 (I-N) DOUBLE PRECISION DTAB(3,81) DOUBLE PRECISION DTAB1(3,20),DTAB2(3,20),DTAB3(3,20),DTAB4(3,20), 1 DTAB5(3,1) EQUIVALENCE (DTAB(1,1), DTAB1(1,1)),(DTAB(1,21),DTAB2(1,1)), 1 (DTAB(1,41),DTAB3(1,1)),(DTAB(1,61),DTAB4(1,1)), 2 (DTAB(1,81),DTAB5(1,1)) SAVE NTAB,IWARN,ITAB,DTAB CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DTAB(1,I) is the year (in decimal form), C C DTAB(2,I) is the JULIAN date (in days), and C C DTAB(3,I) is the difference ET minus UTC, or TDT minus UTC C C in sec, as taken from the Astronomical Ephemeris or the Bulletin C C of the International Earth Rotation Servive. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DATA DTAB1/ 1955.50000D0, 2435289.50000D0, 31.590D0, 1 1956.50000D0, 2435655.50000D0, 32.060D0, 2 1957.50000D0, 2436020.50000D0, 31.820D0, 3 1958.50000D0, 2436385.50000D0, 32.690D0, 4 1959.50000D0, 2436750.50000D0, 33.050D0, 5 1960.50000D0, 2437116.50000D0, 33.160D0, 6 1961.50000D0, 2437481.50000D0, 33.590D0, 7 1962.00000D0, 2437665.50000D0, 34.032D0, 8 1962.50000D0, 2437846.50000D0, 34.235D0, 9 1963.00000D0, 2438030.50000D0, 34.441D0, * 1963.50000D0, 2438211.50000D0, 34.644D0, 1 1964.00000D0, 2438395.50000D0, 34.950D0, 2 1964.50000D0, 2438577.50000D0, 35.286D0, 3 1965.00000D0, 2438761.50000D0, 35.725D0, 4 1965.50000D0, 2438942.50000D0, 36.160D0, 5 1966.00000D0, 2439126.50000D0, 36.498D0, 6 1966.50000D0, 2439307.50000D0, 36.968D0, 7 1967.00000D0, 2439491.50000D0, 37.444D0, 8 1967.50000D0, 2439672.50000D0, 37.913D0, 9 1968.00000D0, 2439856.50000D0, 38.390D0/ DATA DTAB2/ 1968.25000D0, 2439947.50000D0, 38.526D0, 1 1968.50000D0, 2440038.50000D0, 38.760D0, 2 1968.75000D0, 2440130.50000D0, 39.000D0, 3 1969.00000D0, 2440222.50000D0, 39.238D0, 4 1969.25000D0, 2440312.50000D0, 39.472D0, 5 1969.50000D0, 2440403.50000D0, 39.707D0, 6 1969.75000D0, 2440495.50000D0, 39.946D0, 7 1970.00000D0, 2440587.50000D0, 40.185D0, 8 1970.25000D0, 2440677.50000D0, 40.420D0, 9 1970.50000D0, 2440768.50000D0, 40.654D0, * 1970.75000D0, 2440860.50000D0, 40.892D0, 1 1971.00000D0, 2440952.50000D0, 41.131D0, 2 1971.08500D0, 2440983.50000D0, 41.211D0, 3 1971.16200D0, 2441011.50000D0, 41.284D0, 4 1971.24700D0, 2441042.50000D0, 41.364D0, 5 1971.32900D0, 2441072.50000D0, 41.442D0, 6 1971.41400D0, 2441103.50000D0, 41.522D0, 7 1971.49600D0, 2441133.50000D0, 41.600D0, 8 1971.58100D0, 2441164.50000D0, 41.680D0, 9 1971.66600D0, 2441195.50000D0, 41.761D0/ DATA DTAB3/ 1971.74800D0, 2441225.50000D0, 41.838D0, 1 1971.83300D0, 2441256.50000D0, 41.919D0, 2 1971.91500D0, 2441286.50000D0, 41.996D0, 3 1971.99999D0, 2441317.49999D0, 42.184D0, 4 1972.00000D0, 2441317.50000D0, 42.184D0, 5 1972.49999D0, 2441499.49999D0, 42.184D0, 6 1972.50000D0, 2441499.50000D0, 43.184D0, 7 1972.99999D0, 2441683.49999D0, 43.184D0, 8 1973.00000D0, 2441683.50000D0, 44.184D0, 9 1973.99999D0, 2442048.49999D0, 44.184D0, * 1974.00000D0, 2442048.50000D0, 45.184D0, 1 1974.99999D0, 2442413.49999D0, 45.184D0, 2 1975.00000D0, 2442413.50000D0, 46.184D0, 3 1975.99999D0, 2442778.49999D0, 46.184D0, 4 1976.00000D0, 2442778.50000D0, 47.184D0, 5 1976.99999D0, 2443144.49999D0, 47.184D0, 6 1977.00000D0, 2443144.50000D0, 48.184D0, 7 1977.99999D0, 2443509.49999D0, 48.184D0, 8 1978.00000D0, 2443509.50000D0, 49.184D0, 9 1978.99999D0, 2443874.49999D0, 49.184D0/ DATA DTAB4/ 1979.00000D0, 2443874.50000D0, 50.184D0, 1 1979.99999D0, 2444239.49999D0, 50.184D0, 2 1980.00000D0, 2444239.50000D0, 51.184D0, 3 1981.49999D0, 2444786.49999D0, 51.184D0, 4 1981.50000D0, 2444786.50000D0, 52.184D0, 5 1982.49999D0, 2445151.49999D0, 52.184D0, 6 1982.50000D0, 2445151.50000D0, 53.184D0, 7 1983.49999D0, 2445516.49999D0, 53.184D0, 8 1983.50000D0, 2445516.50000D0, 54.184D0, 9 1985.49999D0, 2446247.49999D0, 54.184D0, * 1985.50000D0, 2446247.50000D0, 55.184D0, 1 1987.99999D0, 2447161.49999D0, 55.184D0, 2 1988.00000D0, 2447161.50000D0, 56.184D0, 3 1989.99999D0, 2447892.49999D0, 56.184D0, 4 1990.00000D0, 2447892.50000D0, 57.184D0, 5 1990.99999D0, 2448257.49999D0, 57.184D0, 6 1991.00000D0, 2448257.50000D0, 58.184D0, 7 1992.49999D0, 2448804.49999D0, 58.184D0, 8 1992.50000D0, 2448804.50000D0, 59.184D0, 9 1993.49999D0, 2449169.49999D0, 59.184D0/ DATA DTAB5/ 1993.50000D0, 2449169.50000D0, 60.184D0/ DATA NTAB/81/,IWARN/1/,ITAB/1/ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Print table for the first call of ETMUTC. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(IPRINT.NE.2) GOTO 10 IF(ITAB.EQ.0) GOTO 10 WRITE(IUN6,7001) DO 1 I=1,NTAB WRITE(IUN6,7002) I,DTAB(1,I),DTAB(2,I),DTAB(3,I) 1 CONTINUE ITAB=0 10 CONTINUE IF(DTUJD.LT.DTAB(2,NTAB)) GOTO 11 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DTUJD exceeds last tabulated epoch DTAB(NTAB,3). C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 12 DDT=DTAB(3,NTAB) IF(IWARN.EQ.1) WRITE(IUN6,7051) DTAB(1,NTAB) IF(IWARN.EQ.1) WRITE(*,7051) DTAB(1,NTAB) IWARN=0 RETURN 11 DO 20 I=1,NTAB IF(DTUJD-DTAB(2,I)) 21,22,20 22 DDT=DTAB(3,I) RETURN 21 N=I-1 GOTO 23 20 CONTINUE 23 DDT=(DTAB(3,N+1)*(DTUJD-DTAB(2,N))-DTAB(3,N)*(DTUJD-DTAB(2,N+1)))/ 1 (DTAB(2,N+1)-DTAB(2,N)) RETURN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Format statements: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 7001 FORMAT(//' Routine ETMUTC, version 930712 FTN 77.'// 1' List of tables:'// 2' No. Juld DTX DTY'//) 7002 FORMAT(I10,2F15.5,F10.3) 7051 FORMAT(/' ***** Warning from routine ETMUTC, version 930712.'/ 1' ***** Epoch exceeds the last tabulated value:',F10.5/ 2' ***** DDT of last tabulated epoch is used.'/ 3' ***** Please try to update tabels in routine ETMUC.'/) END C SUBROUTINE ETNLPF(IUN6,KFILT,NFI,DLF,CFILT) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Routine ETNLPF, version 930628 FORTRAN77. C C C C The routine defines symmetrical nor-recursive numerical lowpass C C filters (FIR) to be used within Earth tide analysis program C C ETERNA. C C C C Reference: PERTSEV, B. 1957: On the calculation of drift curve C C in observation of bodily tides. Bulletin C C d' Informations, Marees Terrestres, no. 5, 71-72, C C Bruxelles 1957. C C C C PERTSEV, B. 1959: Ob outchetie spolzaniya nulia C C pir nabloudenij ouprougikh prilivov, Izv. Akad. C C Naouk SSR, no. 4, 1959. C C C C Input parameter description: C C ---------------------------- C C C C IUN6... unit number of formatted line printer file. C C KFILT... Parameter for the selection of a specific numerical C C lowpass filter. C C The default for KFILT is 7 (this seems to be the C C best trade off between length and quality). C C KFILT= 0: no lowpass filter. C C KFILT= 1: lowpass filter of PERTSEV 1957 is used C C with a length of 37 hours. C C KFILT= 2: lowpass filter of PERTSEV 1959 is used C C with a length of 51 hours. C C KFILT= 3: lowpass filter of 143 hours length is C C used, computed by inverse FOURIER- C C transform of the desired frequency C C transfer function of the filter using C C program NUMFIL.FOR at 19880218. C C This filter still exists because of C C continuity but should be replaced by the C C better filter KFILT=7. C C KFILT= 4: lowpass filter of 239 hours length is C C used, computed by inverse FOURIER- C C transform of the desired frequency C C transfer function of the filter using C C program NUMFIL.FOR at 19880218. C C This filter still exists because of C C continuity but should be replaced by the C C better filter KFILT=8. C C KFILT= 5: lowpass filter of 49 hours length is used C C taken from SCHUELLER's HYCON-MC analysis C C program. C C KFILT= 6: lowpass filter of 49 hours length, C C computed by least squares adjustment of C C the desired frequency transfer function C C by WENZEL at 930628 using program C C SYMFIL.FOR, variant 005. C C KFILT=7: bandreject filter of 145 hours length, C C rejection between 0.8 and 5.0 cpd, gives C C as inverse filter a bandpass between 0.8 C C and 5.0 cpd. Computed by least squares C C adjustment of the desired frequency C C transfer function by WENZEL at 930628 C C using program SYMFIL.FOR, variant 007. C C KFILT=8: bandreject filter of 241 hours length, C C rejection between 1.0 and 5.0 cpd, gives C C as inverse filter a bandpass between 1.0 C C and 5.0 cpd. Computed by least squares C C adjustment of the desired frequency C C transfer function by WENZEL at 930627 C C using program SYMFIL.FOR, variant 008. C C C C The quality of the different filters is described by the C C following table: C C C C KFILT length damping [dB] damping [dB] ripple [dB] C C below 0.1 cpd at 0.5 cpd 0.85... 5.0 cpd C C C C 1 37 -30.07 -4.64 2.26 C C 2 51 -63.18 -11.10 4.54 C C 3 143 -56.81 -36.59 1.59 C C 4 239 -69.52 -48.93 0.50 C C 5 49 -64.54 -11.84 1.89 C C 6 49 -69.47 -14.36 1.86 C C 7 145 -153.40 -37.23 0.78 C C 8 241 -167.88 -51.78 0.32 C C C C Output parameter description: C C ----------------------------- C C C C NFI... length of filter in hours (number of filter C C coefficients). C C DLF... array of lowpass filter coefficients (1...239). C C CFILT... name of the filter (CHARACTER*12). C C C C Numerical accuracy: C C ------------------- C C C C The routine has been tested on an IBM-PC using DOUBLE PRECISION C C (i.e. 15 digits) for all non-integer variables. C C C C Routine creation: 880218 by H.-G. Wenzel. C C Geodaetisches Institut, C C Universitaet Karlsruhe, C C Englerstr. 7, C C D-76128 KARLSRUHE 1, C C Germany. C C Tel: 0049-721-6082307, C C FAX: 0049-721-694552. C C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de C C Last modification: 930628 by H.-G. Wenzel. C C**********************************************************************C IMPLICIT DOUBLE PRECISION (D) IMPLICIT INTEGER*4 (I-N) DOUBLE PRECISION DLF(239),DLF1(19),DLF51(26),DLF143(72) DOUBLE PRECISION DL239A(79),DL239B(41),DLF5(25),DLF6(25) DOUBLE PRECISION DLF7(73),DLF8(121) CHARACTER CFILT*12 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Define coefficients of 37 hour numerical lowpass filter from C C PERTSEV 1957, cf. WENZEL 1976, page XXI. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DATA DLF1/ 0.0666666667D0, 0.0000000000D0, 0.0000000000D0, 1 0.0000000000D0, 0.0000000000D0, 0.0666666667D0, 0.0000000000D0, 2 0.0000000000D0, 0.0666666667D0, 0.0000000000D0, 0.0666666667D0, 3 0.0000000000D0, 0.0000000000D0, 0.0666666667D0, 0.0000000000D0, 4 0.0666666667D0, 0.0666666667D0, 0.0000000000D0, 0.0666666667D0/ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Define coefficients of 51 hour numerical lowpass filter from C C PERTSEV 1959, cf. WENZEL 1976, page XXI. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DATA DLF51/ -0.0386100386D0,-0.0386100386D0, 0.0000000000D0, 1 0.0000000000D0, 0.0000000000D0, 0.0000000000D0, 0.0000000000D0, 2 0.0666666667D0, 0.0000000000D0, 0.0000000000D0, 0.0000000000D0, 3 0.0000000000D0, 0.0666666667D0, 0.0000000000D0, 0.0000000000D0, 4 0.0666666667D0, 0.0000000000D0, 0.0666666667D0, 0.0000000000D0, 5 0.0000000000D0, 0.0666666667D0, 0.0000000000D0, 0.0666666667D0, 6 0.0666666667D0, 0.0386100386D0, 0.1438867439D0/ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Define coefficients of 143 hour numerical lowpass filter. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DATA DLF143/ 0.0003157695D0, 0.0002850986D0, 0.0002439589D0, 1 0.0001913863D0, 0.0001264130D0, 0.0000482874D0,-0.0000432812D0, 2-0.0001478733D0,-0.0002641261D0,-0.0003895397D0,-0.0005203443D0, 3-0.0006514461D0,-0.0007764701D0,-0.0008879070D0,-0.0009773704D0, 4-0.0010359605D0,-0.0010547247D0,-0.0010251972D0,-0.0009399940D0, 5-0.0007934331D0,-0.0005821445D0,-0.0003056342D0, 0.0000332385D0, 6 0.0004279060D0, 0.0008678952D0, 0.0013388165D0, 0.0018225602D0, 7 0.0022977131D0, 0.0027401971D0, 0.0031241197D0, 0.0034228162D0, 8 0.0036100510D0, 0.0036613339D0, 0.0035553004D0, 0.0032750948D0, 9 0.0028096936D0, 0.0021551009D0, 0.0013153517D0, 0.0003032613D0, *-0.0008591336D0,-0.0021404843D0,-0.0035004235D0,-0.0048902555D0, 1-0.0062540453D0,-0.0075300873D0,-0.0086527186D0,-0.0095544247D0, 2-0.0101681712D0,-0.0104298800D0,-0.0102809599D0,-0.0096707954D0, 3-0.0085590926D0,-0.0069179844D0,-0.0047338019D0,-0.0020084284D0, 4 0.0012398342D0, 0.0049759384D0, 0.0091483392D0, 0.0136898276D0, 5 0.0185189400D0, 0.0235419109D0, 0.0286551118D0, 0.0337479045D0, 6 0.0387058123D0, 0.0434139046D0, 0.0477602767D0, 0.0516395027D0, 7 0.0549559358D0, 0.0576267397D0, 0.0595845368D0, 0.0607795779D0, 8 0.0611813522D0/ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Define coefficients of 239 hour numerical lowpass filter. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DATA DL239A/ -0.0000003459D0,-0.0000011382D0,-0.0000019040D0, 1-0.0000020720D0,-0.0000010321D0, 0.0000017969D0, 0.0000068966D0, 2 0.0000145807D0, 0.0000249356D0, 0.0000377724D0, 0.0000525981D0, 3 0.0000686066D0, 0.0000846941D0, 0.0000994991D0, 0.0001114664D0, 4 0.0001189331D0, 0.0001202331D0, 0.0001138150D0, 0.0000983666D0, 5 0.0000729394D0, 0.0000370651D0,-0.0000091452D0,-0.0000649241D0, 6-0.0001288027D0,-0.0001986023D0,-0.0002714657D0,-0.0003439321D0, 7-0.0004120544D0,-0.0004715573D0,-0.0005180302D0,-0.0005471490D0, 8-0.0005549170D0,-0.0005379139D0,-0.0004935403D0,-0.0004202447D0, 9-0.0003177209D0,-0.0001870601D0,-0.0000308496D0, 0.0001467947D0, * 0.0003402671D0, 0.0005425882D0, 0.0007455980D0, 0.0009402242D0, 1 0.0011168179D0, 0.0012655482D0, 0.0013768396D0, 0.0014418373D0, 2 0.0014528790D0, 0.0014039540D0, 0.0012911254D0, 0.0011128947D0, 3 0.0008704869D0, 0.0005680375D0, 0.0002126658D0,-0.0001855784D0, 4-0.0006139007D0,-0.0010570830D0,-0.0014979349D0,-0.0019178663D0, 5-0.0022975640D0,-0.0026177523D0,-0.0028600092D0,-0.0030076097D0, 6-0.0030463602D0,-0.0029653905D0,-0.0027578660D0,-0.0024215850D0, 7-0.0019594302D0,-0.0013796444D0,-0.0006959079D0, 0.0000727977D0, 8 0.0009025499D0, 0.0017650676D0, 0.0026285091D0, 0.0034584659D0, 9 0.0042191268D0, 0.0048745789D0, 0.0053902005D0, 0.0057341020D0/ DATA DL239B/ 0.0058785602D0, 0.0058013927D0, 0.0054872170D0, 1 0.0049285412D0, 0.0041266363D0, 0.0030921468D0, 0.0018454020D0, 2 0.0004164023D0,-0.0011555371D0,-0.0028224882D0,-0.0045289516D0, 3-0.0062132168D0,-0.0078090073D0,-0.0092473692D0,-0.0104587464D0, 4-0.0113751813D0,-0.0119325702D0,-0.0120729000D0,-0.0117463900D0, 5-0.0109134640D0,-0.0095464807D0,-0.0076311569D0,-0.0051676267D0, 6-0.0021710905D0, 0.0013279801D0, 0.0052841008D0, 0.0096375005D0, 7 0.0143153041D0, 0.0192331451D0, 0.0242971663D0, 0.0294063546D0, 8 0.0344551412D0, 0.0393361916D0, 0.0439433008D0, 0.0481743039D0, 9 0.0519339131D0, 0.0551363910D0, 0.0577079759D0, 0.0595889822D0, * 0.0607355080D0, 0.0611206964D0/ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Define coefficients of 49 hour numerical lowpass filter, C C taken from SCHUELLER's HYCON-MC analysis program (this filter is C C not very good, it is implemented just for comparison purpose). C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DATA DLF5/ -0.0276110400D0,-0.0216086400D0,-0.0158616600D0, 1-0.0103701100D0,-0.0051339680D0,-0.0001532528D0, 0.0045720420D0, 2 0.0090419150D0, 0.0132563700D0, 0.0172154000D0, 0.0209190100D0, 3 0.0243671900D0, 0.0275599600D0, 0.0304973100D0, 0.0331792300D0, 4 0.0356057300D0, 0.0377768100D0, 0.0396924700D0, 0.0413527100D0, 5 0.0427575300D0, 0.0439069200D0, 0.0448009000D0, 0.0454394500D0, 6 0.0458225800D0, 0.0459503000D0/ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Define coefficients of 49 hour numerical lowpass filter, C C computed by least squares adjustment of the desired frequency C C transfer function of the filter using program SYMFIL.FOR, C C at 930628, variant 005. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DATA DLF6/ -.014137780885889D0, -.013624837420017D0, 1 -.011216507417872D0, -.009898073147066D0, -.006978605732010D0, 2 -.004795964664175D0, -.001405243674584D0, .001599349619534D0, 3 .005379061497206D0, .009062534703677D0, .013105895320699D0, 4 .017238380257761D0, .021379527407788D0, .025668653877614D0, 5 .029707398051727D0, .033830864233968D0, .037543679061798D0, 6 .041184461877835D0, .044340166115800D0, .047220017922438D0, 7 .049599939018570D0, .051506418937014D0, .052926652116268D0, 8 .053731562876735D0, .054064900090366D0/ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Define coefficients of 145 hour numerical bandreject filter, C C computed by least squares adjustment of the desired frequency C C transfer function of the filter using program SYMFIL.FOR, C C at 930628, variant 007. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DATA (DLF7(J),J=1,57) / 1 -.001382847852496D0, -.002978380523749D0, -.000920253794430D0, 2 .002425751898133D0, .002186129435109D0, .000448584004330D0, 3 .001501281378205D0, .003663628113565D0, .003024503608685D0, 4 .000414992641781D0, -.000528184156354D0, .000726939607111D0, 5 .000920078864678D0, -.001385330263747D0, -.003382385854319D0, 6 -.002892033875532D0, -.001836609200941D0, -.002566782910910D0, 7 -.004043931644741D0, -.004027103855561D0, -.002674782062247D0, 8 -.001814319802053D0, -.001787934255117D0, -.001419201265932D0, 9 -.000351498940939D0, .000735343689753D0, .001779772083634D0, * .003080725054110D0, .004008776400498D0, .003960528252744D0, 1 .004116839769344D0, .005793965333947D0, .007372659808072D0, 2 .006198096678327D0, .003539825675513D0, .003550476248496D0, 3 .006144188736919D0, .005864601289585D0, .000627958784182D0, 4 -.003444614332327D0, -.001220833508571D0, .001895625270243D0, 5 -.002469464648339D0, -.011048409879104D0, -.012326591760107D0, 6 -.005409997644967D0, -.003473822574344D0, -.012859882534529D0, 7 -.020528266826841D0, -.013571006387825D0, -.001800445803859D0, 8 -.004312336812885D0, -.017580343137789D0, -.017782074467984D0, 9 .001068549953901D0, .013605424960193D0, .001923541512108D0/ DATA (DLF7(J),J=58,73) / 1 -.011882371338061D0, .002830618378610D0, .033151557400089D0, 2 .036635625824296D0, .009721344108865D0, .002191069063451D0, 3 .041199093952236D0, .077335262788393D0, .052085916612365D0, 4 -.000176378010294D0, .016680188941672D0, .110158097881496D0, 5 .144434185321490D0, -.002252796303815D0, -.252236167435206D0, 6 .622731268679572D0/ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Define coefficients of 241 hour numerical bandreject filter, C C computed by least squares adjustment of the desired frequency C C transfer function of the filter using program SYMFIL.FOR, C C at 930628, variant 008. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DATA (DLF8(J),J=1,57)/ 1 .000290267264581D0, .000361404738067D0, .000058498741689D0, 2 -.000006593633207D0, -.000184830919523D0, -.000623441163443D0, 3 -.000723760167717D0, -.000276448608695D0, -.000064710571002D0, 4 -.000462474665584D0, -.000595047927989D0, .000106501571412D0, 5 .000738584451467D0, .000482499708227D0, .000113888923717D0, 6 .000591522704785D0, .001262360599822D0, .001016287195346D0, 7 .000281074195059D0, .000215479845356D0, .000655853123205D0, 8 .000547049740554D0, -.000114777572607D0, -.000409642396305D0, 9 -.000209436314139D0, -.000237600082115D0, -.000647919253087D0, * -.000928113563303D0, -.001022945538863D0, -.001282319161605D0, 1 -.001503740911559D0, -.001233742772022D0, -.000785587586708D0, 2 -.000780606025348D0, -.000947803264414D0, -.000572249800445D0, 3 .000139823670465D0, .000386238991944D0, .000232824574474D0, 4 .000471591058126D0, .001134278696183D0, .001429142785714D0, 5 .001197329127838D0, .001226954270461D0, .001902028630981D0, 6 .002545912347880D0, .002430955891520D0, .001729036343091D0, 7 .001157816386045D0, .001095157789585D0, .001216791216624D0, 8 .000893177118717D0, .000017746595035D0, -.000703459058713D0, 9 -.000743455318216D0, -.000749492255399D0, -.001603759466679D0/ DATA (DLF8(J),J=58,114)/ 1 -.002742111857542D0, -.002847821685452D0, -.002188647618656D0, 2 -.002406401826882D0, -.003654408908016D0, -.004096189416041D0, 3 -.003008879172286D0, -.002094710110113D0, -.002512197869440D0, 4 -.002882543584930D0, -.001815391805948D0, -.000273129719995D0, 5 .000367241464943D0, .000629674872347D0, .001493326699603D0, 6 .002531739533750D0, .003234049802549D0, .004187154339736D0, 7 .005337870344607D0, .005409512046386D0, .004615233135387D0, 8 .005074984867706D0, .006793167992483D0, .006780972597671D0, 9 .004376399793946D0, .003170471000973D0, .004643649178669D0, * .004732146449398D0, .000839259429395D0, -.002837211417737D0, 1 -.002310356969918D0, -.001106548630836D0, -.004065200216209D0, 2 -.008527885305902D0, -.009232975359170D0, -.007819238548454D0, 3 -.009187152198087D0, -.012133819682557D0, -.012345472112452D0, 4 -.011048159434414D0, -.011561647314243D0, -.011459532283642D0, 5 -.007997604120371D0, -.005849510848135D0, -.007994639547899D0, 6 -.006296457513721D0, .004191211764374D0, .010957343122874D0, 7 .004212923389456D0, .000670178884434D0, .018569747820049D0, 8 .038424790264353D0, .028789386560494D0, .007535287259817D0, 9 .023409211872283D0, .066773877885917D0, .067026848814559D0/ DATA (DLF8(J),J=115,121)/ 1 .013871958440154D0, .004172698405786D0, .091317960780299D0, 2 .150333399502551D0, .015689896712377D0, -.256620703639583D0, 3 .600477706780038D0/ IF(KFILT.EQ.0) GOTO 1000 GOTO (1100,1200,1300,1400,1500,1600,1700,1800),KFILT GOTO 1700 1000 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C KFILT = 0, no lowpass filter. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC NFI=1 DLF(1)=0.D0 CFILT='NO FILTER ' RETURN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C KFILT = 1, lowpass filter of PERTSEV 1957, 37 hours length. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1100 CFILT='PERTSEV 1957' NFI=37 NFI2=NFI/2 NFI3=NFI2+1 DO 1110 IF=1,NFI3 DLF(IF)=DLF1(IF) JF=NFI-IF+1 DLF(JF)=DLF1(IF) 1110 CONTINUE GOTO 1900 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C KFILT = 2, lowpass filter of PERTSEV 1959, 51 hours length. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1200 CFILT='PERTZEV 1959' NFI=51 NFI2=NFI/2 NFI3=NFI2+1 DO 1210 IF=1,NFI3 DLF(IF)=DLF51(IF) JF=NFI-IF+1 DLF(JF)=DLF51(IF) 1210 CONTINUE GOTO 1900 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C KFILT = 3, lowpass filter computed by WENZEL 1988, 143 hours. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1300 CFILT='WENZEL 143 ' NFI=143 NFI2=NFI/2 NFI3=NFI2+1 DO 1310 IF=1,NFI3 DLF(IF)=DLF143(IF) JF=NFI-IF+1 DLF(JF)=DLF143(IF) 1310 CONTINUE GOTO 1900 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C KFILT = 4, lowpass filter computed by WENZEL 1988, 239 hours. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1400 CFILT='WENZEL 239 ' NFI=239 NFI2=NFI/2 NFI3=NFI2+1 DO 1410 IF=1,79 DLF(IF)=DL239A(IF) JF=NFI-IF+1 DLF(JF)=DL239A(IF) 1410 CONTINUE DO 1420 IF=1,41 JF=IF+79 DLF(JF)=DL239B(IF) JF=NFI-JF+1 DLF(JF)=DL239B(IF) 1420 CONTINUE GOTO 1900 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C KFILT = 5, define numerical lowpass filter taken from program C C HYCON-MC, length of 49 hours. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1500 CFILT='SCHUELLER 78' NFI=49 NFI2=NFI/2 NFI3=NFI2+1 DO 1510 IF=1,NFI3 DLF(IF)=DLF5(IF) JF=NFI-IF+1 DLF(JF)=DLF5(IF) 1510 CONTINUE GOTO 1900 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C KFILT = 6, define numerical lowpass filter with 49 hours length, C C computed by least squares adjsutment of the desired C C frequency transfer function of the filter using C C SYMFIL.FOR, at 930628, variant 005. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1600 CFILT='WENZEL 49 ' NFI=49 NFI2=NFI/2 NFI3=NFI2+1 DO 1610 IF=1,NFI3 DLF(IF)=DLF6(IF) JF=NFI-IF+1 DLF(JF)=DLF6(IF) 1610 CONTINUE GOTO 1900 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C KFILT = 7, define numerical bandreject filter with 145 hours C C length, computed by least squares adjsutment of the C C desired frequency transfer function of the filter C C using program SYMFIL.FOR, at 930627, variant 008. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1700 CFILT='WENZEL 145' NFI=145 NFI2=NFI/2 NFI3=NFI2+1 DO 1710 IF=1,NFI3 DLF(IF)=DLF7(IF) JF=NFI-IF+1 DLF(JF)=DLF7(IF) 1710 CONTINUE GOTO 1900 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C KFILT = 8, define numerical bandreject filter with 241 hours C C length, computed by least squares adjsutment of the C C desired frequency transfer function of the filter C C using program SYMFIL.FOR, at 930627, variant 008. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1800 CFILT='WENZEL 241' NFI=241 NFI2=NFI/2 NFI3=NFI2+1 DO 1810 IF=1,NFI3 DLF(IF)=DLF8(IF) JF=NFI-IF+1 DLF(JF)=DLF8(IF) 1810 CONTINUE 1900 CONTINUE WRITE(IUN6,7001) CFILT,NFI WRITE(IUN6,7006) RETURN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Format statements: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 7001 FORMAT(' Routine ETNLPF, version 930628 FTN77.'/ 1 ' Used numerical filter is ',A12/ 2 ' Length of the used numerical filter is ',I6,' hours') 7006 FORMAT(' Routine ETNLPF finished the execution.'/) END C SUBROUTINE ETPOTA(IUN4,IUN6,IUN14,IPRINT,IMODEL,DLAT,DLON,DH, 1 DGRAV,DAZ,IC,ITY,ITM,ITD,DTH,DDT0,MAXNW,NRW,DTHAM,DTHPH,DTHFR, 2 DBODY,NW) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Routine ETPOTA, version 930710 FORTRAN 77. C C C C The routine ETPOTA computes amplitudes, phases, frequencies and C C body tide amplitude factors for a number of different Earth tide C C components using three different tidal potential developments. C C C C Attention: This routine has finally not been tested for vertical C C and horizontal displacements and for shear tidal C C strain !!!! C C C C All variables with D as first character are DOUBLE PRECISION. C C C C Input parameter description: C C ---------------------------- C C C C IUN4... formatted unit, on which the tidal potential C C development has to be stored before the execution C C of routine ETPOTA (e.g. file ETCPOT.DAT). C C IUN6... formatted line printer unit. C C IUN14... unformatted copy of IUN4. This unit will be opened C C as file ETCPOT.UFT during the execution of routine C C ETPOTA with STATUS=OLD if it exists and with C C STATUS=NEW, if it does not exist. If the STATUS is C C NEW, ETCPOT.DAT will be established during the C C execution of routine ETPOTA. C C IPRINT... printout parameter. C C for IPRINT = 0, nothing will be printed. C C for IPRINT = 1, a short list will be printed. C C for IPRINT = 2, a long list will be printed C C (including the tidal potential development). C C IMODEL... parameter for selecting the tidal potential C C development. C C IMODEL = 0: DOODSON 1921 tidal potential develop- C C ment with 378 waves. C C IMODEL = 1: CARTWRIGHT-TAYLOR-EDDEN 1973 tidal C C potential development with 505 waves. C C IMODEL = 2: TAMURA 1987 tidal potential develop- C C ment with 1200 waves. C C IMODEL = 3: BUELLESFELD 1985 tidal potential C C development with 656 waves. C C DLAT... ellipsoidal latitude referring to Geodetic C C Reference System 1980 in degree. C C DLON... ellipsoidal longitude referring to Geodetic C C Reference System 1980 in degree, positive east of C C Greenwhich. C C DH... ellipsoidal height referring to Geodetic Reference C C System 1980 in meter. C C DGRAV... gravity in m/s**2. If the gravity is input below C C 1 m/s**2, the gravity will be replaced by the C C computed normal gravity for reference system GRS80. C C DAZ... azimuth in degree from north direction (only valid C C for tidal tilt, horizontal displacement, and C C horizontal strain). C C IC... Earth tide component to be computed. C C IC=-1: tidal potential in m**2/s**2. C C IC= 0: vertical tidal acceleration (gravity tide), C C in nm/s**2 (positive downwards). C C IC= 1: horizontal tidal acceleration (tidal tilt) C C in azimuth DAZ in mas = arc sec/1000. C C IC= 2: vertical tidal displacement, geodetic C C coefficients in mm (positive upwards). C C IC= 3: horizontal tidal displacement in azimuth C C DAZ in mm. C C IC= 4: vertical tidal strain in 10**-9 = nstr. C C IC= 5: horizontal tidal strain in azimuth DAZ C C in 10**-9 = nstr. C C IC= 6: areal tidal strain in 10**-9 = nstr. C C IC= 7: shear tidal strain in 10**-9 = nstr. C C IC= 8: volume tidal strain in 10**-9 = nstr. C C IC= 9: ocean tides, geodetic coefficients in C C millimeter. C C IR... printout parameter for the tidal potential development. C C IR= 0: no printout of the tidal potential development. C C IR= 1: tidal potential development and the development C C of the specific tidal component will be printed. C C ITY... year of initial epoch of tidal development. C C ITM... month of initial epoch of tidal development. C C ITD... day of initial epoch of tidal development. C C DTH... hour of initial epoch of tidal development (UTC). C C Example: September 27th 1984, 10 o'clock UTC is C C ITY = 1984, ITM = 9, ITD = 27, DTH = 10.D0 C C MAXNW... maximum number of waves, used for DIMENSION of C C arrays NRW, DTHAM, DTHPH, DTHFR, DBODY. Although C C the TAMURA 1987 potential development contains 1200 C C waves only, the CTED 1973 and the BUELLESFELD 1985 C C tidal potential developments contain some few waves C C with small amplitude, which do not exist in the C C TAMURA 1987 potentail development. Thus, the total C C maximum number of waves MAXNW of file ETCPOT.DAT is C C 1214. C C C C Output parameter description: C C ----------------------------- C C C C DDT0... difference TDT - UTC at initial epoch in sec. C C NRW... array(1...MAXNW) of wave numbers. C C DTHAM... array (1...MAXNW) of tidal amplitudes given in C C units of CUNIT(IC2). C C DTHPH... array (1...MAXNW) of tidal phases in radians at C C initial epoch. C C DTHFR... array (1...MAXNW) of tidal frequencies in radian C C per hour. C C DBODY... array (1...MAXNW) of bady tide amplitude factors C C for tidal gravity and tidal tilt. In order to C C compute the body tide, the amplitudes DTHAM have to C C be multiplied by DBODY. In case of IRIGID=1, all C C body tide amplitude factors are set to 1.0000 C C NW... number of defined tidal waves. C C C C Used routines: C C -------------- C C C C ETASTE: computes astronomical elements. C C ETGCOF: computes geodetic coefficients. C C ETJULD: computes Julian date. C C ETLOVE: computes latitude dependent elastic parameters (called C C ETGCOF). C C ETMUTC: computes the difference TDT minus UTC (called by ETASTE).C C GEOEXT: computes jobtime. C C C C Numerical accuracy: C C ------------------- C C C C The routine has been tested under operation systems UNIX and C C MS-DOS with 15 digits in DOUBLE PRECISION. C C C C Routine creation: 880427 by H.-G. Wenzel. C C Geodaetisches Institut, C C Universitaet Karlsruhe, C C Englerstr. 7, C C D-76128 KARLSRUHE 1, C C Germany. C C Tel: 0049-721-6082307, C C FAX: 0049-721-694552. C C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de C C Last modification: 930706 by H.-G. Wenzel. C C**********************************************************************C IMPLICIT DOUBLE PRECISION (D) IMPLICIT INTEGER*4 (I-N) LOGICAL*1 LEX14 CHARACTER CHEAD(8)*10,CENDH*10,CUNIT(11)*8,CWN*4,CMODEL(4)*13 INTEGER*2 NS(8) DOUBLE PRECISION DX(3),DHH(4) DOUBLE PRECISION DAS(8),DASP(8),DGK(12),DPK(12) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C The following DIMENSION statement is concerning the number of C C waves of the tidal potential development, which is MAXNW. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC INTEGER*4 NRW(MAXNW) DOUBLE PRECISION DTHAM(MAXNW),DTHPH(MAXNW),DTHFR(MAXNW), 1 DBODY(MAXNW) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C The following DIMENSION statement is concerning the elastic C C Earth model for the different degree and order constituents. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DOUBLE PRECISION DELTA(12) COMMON /UNITS/ CUNIT,IC2 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C COMMON /CONST/: C C DPI... 3.1415.... C C DPI2... 2.D0*DPI C C DRAD... DPI/180.D0 C C DRO... 180.D0/DPI C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC COMMON /CONST/ DPI,DPI2,DRAD,DRO CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C COMMON /LOVE/ contains gravimeter factors, LOVE-numbers, SHIDA- C C numbers and tilt factors for degree 2...4 at latitude DLAT: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DIMENSION DGLAT(12),DHLAT(12),DKLAT(12),DLLAT(12),DTLAT(12) COMMON /LOVE/ DOM0,DOMR,DGLAT,DGR,DHLAT,DHR,DKLAT,DKR,DLLAT,DLR, 1 DTLAT,DTR DATA CENDH/'C*********'/ DATA CMODEL/'DOODSON 1921 ','CTED 1973 ','TAMURA 1987 ', 1 'BUELLESF.1985'/ IF(IPRINT.EQ.0) GOTO 10 WRITE(IUN6,7001) CMODEL(IMODEL+1) WRITE(*,7001) CMODEL(IMODEL+1) 10 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Test, whether there exist already unformatted file ETCPOT.UFT: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC OPEN(UNIT=IUN14,FILE='ETCPOT.UFT',FORM='UNFORMATTED',STATUS='OLD', 1 ERR=11) LEX14=.TRUE. REWIND IUN14 GOTO 12 11 OPEN(UNIT=IUN14,FILE='ETCPOT.UFT',FORM='UNFORMATTED',STATUS='NEW') LEX14=.FALSE. REWIND IUN4 12 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute geodetic coefficients and body tide amplitude factors C C for the WAHR-DEHANT-ZSCHAU model. The NDFW resonance is C C approximated by C C C C G0 - GR*(DOM - DOM0)/(DOMR - DOM), C C C C similar equations hold for the other components. C C C C Gravimetric amplitude factors, LOVE numbers h and k for zero to C C third degree tidal potential have been taken from DEHANT 1987, C C table 7, 8 and 9 for elliptical, uniformly rotating, oceanless C C Earth with liquid outer core and inelastic mantle (PREM Earth C C model with inelastic mantle from ZSCHAU) and for the fourth C C degree from DEHANT et al. 1989, table 6). The resonance factors C C GR have been computed to fit the difference between body tide C C amplitude factors at waves O1 and PSI1 from DEHANT 1987, PREM C C model with elastic mantle (table 1...3). The NDFW resonance C C frequency is 15.073729 degree per hour = 1.004915267 CPD UT, C C taken from WAHR 1981 (because it is not given in any of DEHANT's C C papers). C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL ETGCOF(IUN6,IPRINT,DLAT,DLON,DH,DGRAV,DAZ,IC,DGK,DPK) IC2=IC+2 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Define default body tide amplitude factors for components C C IC=2...9. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 50 I=1,12 50 DELTA(I)=1.D0 DELTAR=0.D0 GOTO (100,200,300),IC2 GOTO 1000 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IC=-1, compute body tide amplitude factors for tidal potential: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 100 CONTINUE DO 110 I=1,12 110 DELTA(I)=DKLAT(I) DELTAR=DKR GOTO 1000 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IC=0, compute body tide amplitude factors for vertical component C C (gravity tides): C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 200 CONTINUE DO 210 I=1,12 210 DELTA(I)=DGLAT(I) DELTAR=DGR GOTO 1000 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IC=1: compute body tide amplitude factors for horizontal C C component (tidal tilt): C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 300 CONTINUE DO 310 I=1,12 310 DELTA(I)=DTLAT(I) DELTAR=DKR-DHR 1000 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C compute JULIAN date for initial epoch ITY,ITM,ITD,DTH: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL ETJULD(IUN6,ITY,ITM,ITD,DTH,DT) DT2000=(DT-2451544.D0)/36525.0D0 DT=(DT-2415020.0D0)/36525.0D0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute astronomical elements for initial epoch: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL ETASTE(IUN6,IPRINT,IMODEL,DLON,ITY,ITM,ITD,DTH,DAS,DASP, 1 DDT0) IC2=IC+2 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Interpolation factors for CTED 1973 potential coefficients at C C start epoch by linear least squares interpolation (see WENZEL C C 1976). C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DX(1)=0.550376D0-1.173312D0*DT DX(2)=0.306592D0+0.144564D0*DT DX(3)=0.143033D0+1.028749D0*DT CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Read file header of tidal potential file on unit IUN4: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(LEX14) READ(IUN14) (CHEAD(I),I=1,8) IF(.NOT.LEX14) READ(IUN4,7028) (CHEAD(I),I=1,8) IF(.NOT.LEX14) WRITE(IUN14) (CHEAD(I),I=1,8) WRITE(IUN6,7029) (CHEAD(I),I=1,8) WRITE(*,7029) (CHEAD(I),I=1,8) 1100 IF(LEX14) READ(IUN14) (CHEAD(I),I=1,8) IF(.NOT.LEX14) READ(IUN4,7028) (CHEAD(I),I=1,8) IF(.NOT.LEX14) WRITE(IUN14) (CHEAD(I),I=1,8) IF(IPRINT.GT.1) WRITE(IUN6,7029) (CHEAD(I),I=1,8) IF(IPRINT.GT.1) WRITE(*,7029) (CHEAD(I),I=1,8) IF(CHEAD(1).NE.CENDH) GOTO 1100 IF(LEX14) READ(IUN14) NW IF(.NOT.LEX14) READ(IUN4,7005) NW IF(.NOT.LEX14) WRITE(IUN14) NW IF(NW.GT.MAXNW) GOTO 5000 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Compute tidal development for the specific component from tidal C C potential development: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DSH1=0.D0 DSH2=0.D0 DSH3=0.D0 DSHD=0.D0 DSHT=0.D0 IW=1 NAMPL=0 1110 CONTINUE 1120 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C NRT is the wave number of TAMURA's 1987 tidal potential C C development, NRC is the wave number of CARTWRIGHT-TAYLER-EDDEN C C 1973 tidal potential development. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(LEX14) READ(IUN14,END=2000) NRI,K,(NS(J),J=1,8),NP, 1 (DHH(J),J=1,3),CWN,DHD,DHT,DHTD,DHB IF(.NOT.LEX14) READ(IUN4,7006,END=2000) NRI,K,(NS(J),J=1,8),NP, 1 (DHH(J),J=1,3),CWN,DHD,DHT,DHTD,DHB IF(.NOT.LEX14) WRITE(IUN14) NRI,K,(NS(J),J=1,8),NP, 1 (DHH(J),J=1,3),CWN,DHD,DHT,DHTD,DHB DSHT=DSHT+DHT J=IW-1 C IF(MOD(J,51).EQ.0.AND.IPRINT.EQ.2) WRITE(IUN6,7007) ITY,ITM,ITD, C 1 DTH,IC,NW,CUNIT(IC2) DHH(4)=0.D0 DSH1=DSH1+DHH(1) DSH2=DSH2+DHH(2) DSH3=DSH3+DHH(3) DSHD=DSHD+DHD CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CTED 1973 tidal potential development is used: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 1130 J=1,3 1130 DHH(4)=DHH(4)+DHH(J)*DX(J) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DOODSON's 1921 tidal potential development is used: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(IMODEL.EQ.0) DHH(4)=DHD CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C TAMURA's 1987 tidal potential development is used: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(IMODEL.EQ.2) DHH(4)=DHT+DHTD*DT2000 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C BUELLESFELD's 1985 tidal potential development is used: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(IMODEL.EQ.3) DHH(4)=DHB IF(DABS(DHH(4)).LT.1.D-6) GOTO 1120 NAMPL=NAMPL+1 DC2=0.D0 DC3=0.D0 DO 1140 J=1,8 DC2=DC2+NS(J)*DAS(J) 1140 DC3=DC3+NS(J)*DASP(J) J=(K+1)*K/2-2+NS(1) DC2=DC2+DPK(J)+NP*90.D0 DC1=DHH(4)*DGK(J) DBODY(IW)=DELTA(J) IF(J.EQ.2) DBODY(IW)=DELTA(J)+DELTAR*(DC3-DOM0)/(DOMR-DC3) IF(DC1.GE.0.D0) GOTO 1160 DC1=-DC1 DC2=DC2-180.D0 1160 DC2=DMOD(DC2,360.D0) IF(DC2.GE.0.D0) GOTO 1170 DC2=DC2+360.D0 GOTO 1160 1170 CONTINUE DTHAM(IW)=DC1 DTHPH(IW)=DC2*DRAD DTHFR(IW)=DC3*DRAD NRW(IW)=NRI C IF(IPRINT.EQ.2) WRITE(IUN6,7008) NRI,K,(NS(J),J=1,8),NP, C 1 (DHH(J),J=1,4),DC1,DC2,DC3,CWN,DHD,DHT,DHTD IF(IPRINT.EQ.2) WRITE(IUN6,7011) DC1,DC2,DC3,CWN,DBODY(IW) IF(IPRINT.EQ.2) WRITE(*,7011) DC1,DC2,DC3,CWN,DBODY(IW) IW=IW+1 GOTO 1110 2000 CONTINUE NW=IW-1 IF(IPRINT.EQ.0) RETURN WRITE(IUN6,7009) DSH1,DSH2,DSH3,DSHD,DSHT WRITE(*,7009) DSH1,DSH2,DSH3,DSHD,DSHT WRITE(IUN6,7010) NW,NAMPL WRITE(*,7010) NW,NAMPL WRITE(IUN6,7030) WRITE(*,7030) RETURN 5000 CONTINUE WRITE(IUN6,7050) NW,MAXNW WRITE(*,7050) NW,MAXNW STOP CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Format statements: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 7001 FORMAT(' Routine ETPOTA, version 930710 FORTRAN 77.'// 1' Tidal component development from tidal potential development.'// 2 1X,A13,' tidal potential development is used.'/) 7005 FORMAT(2I5) 7006 FORMAT(I4,9X,10I2,3F7.5,1X,A4,F7.5,F8.6,F9.6,F8.6) 7007 FORMAT(' Routine ETPOTA, version 930710 FORTRAN 77'// 1' Tidal potential and tidal component development'/ 2' Initial epoch',I6,2(1H.,I2),F5.2/ 3' Component: ',I4/' Number of waves: ',I4// 3' no. argum. argum.nrs. NP H1870 H1924', 4' H1960 H(T) ampl. phase frequency'/ 5 55X,A6,' [deg] [deg/h]'/) 7008 FORMAT(1X,I4,10I2,4F7.5,F8.4,F9.4,F12.8,1X,A4/F7.5,F8.6,F9.6) 7009 FORMAT(//' Sum of potential coefficients :'/5F10.6) 7010 FORMAT(//' Number of waves is :',I6/ 1 ' Number of amplitudes > 0 is:',I6/ ) 7011 FORMAT(3F10.5,2X,A6,2X,F10.6) 7028 FORMAT(8A10) 7029 FORMAT(1X,8A10) 7030 FORMAT(///' ***** Routine ETPOTA finished execution.'/) 7050 FORMAT(/' ***** Error in routine ETPOTA, version 930710.'/ 1' ***** The current number of waves:',I5,' exceeds the maximum', 2' number of waves:',I5/ 3' ***** Routine ETPOTA stops the execution.'/) END C SUBROUTINE ETSDER(IUN6,JB,DAL,IDIML,NO,DATLIM,NDLB) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C ROUTINE ETSDER, VERSION 910408 FORTRAN 77. C C C C === VERSION FOR IBM-AT ======================================== C C C C THE ROUTINE ETSDER SEARCHES FOR DATA ERRORS IN THE EARTH TIDE C C OBSERVATIONS OF THE CURRENT BLOCK. THE DATA TEST FILTER DTF IS C C DESCRIBED IN WENZEL, H.-G. 1976...ZUR GENAUIGKEIT VON GRAVIME- C C TRISCHEN ERDEGEZEITENBEOBACHTUNGEN. WISSENSCHAFT- C C LICHE ARBEITEN DER LEHRSTUEHLE FUER GEODAESIE, C C PHOTOGRAMMETRIE UND KARTOGRAPHIE AN DER TECHNI- C C SCHEN UNIVERSITAET HANNOVER NR. 67, HANNOVER 1976. C C C C ALL VARIABLES WITH D AS FIRST CHARACTER ARE DOUBLE PRECISION. C C C C INPUT PARAMETER DESCRIPTION... C C ============================== C C C C IUN6... FORMATTED LINE PRINTER UNIT. C C JB... NUMBER OF CURRENT BLOCK. C C DAL... ARRAY OF EARTH TIDE OBSERVATIONS (1...IDIML). C C IDIML... DIMENSION OF ARRAY DAL. C C NO... NUMBER OF EARTH TIDE OBSERVATIONS IN THE CURRENT C C BLOCK. C C DATLIM... DATA ERROR LIMIT IN UNITS OF THE EARTH TIDE C C OBSERVATIONS. C C C C OUTPUT PARAMETER DESCRIPTION... C C =============================== C C C C NDLB... NUMBER OF DATA ERRORS SUSPECTED IN THE CURRENT C C BLOCK. C C C C USED ROUTINES... C C ================ C C C C ETGREI: computes GREGORIAN date. C C ETNFIL: does numerical filtering. C C C C Routine creation: 901230 by H.-G. Wenzel. C C Geodaetisches Institut, C C Universitaet Karlsruhe, C C Englerstr. 7, C C D-76128 KARLSRUHE 1, C C Germany. C C Tel: 0049-721-6082307, C C FAX: 0049-721-694552. C C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de C C Last modification: 910408 by H.-G. Wenzel. C C**********************************************************************C IMPLICIT DOUBLE PRECISION (D) IMPLICIT INTEGER*4 (I-N) SAVE DTF INTEGER*4 IRECA(300),IRECE(300),IDATA(300),ITIMA(300),IDATE(300), 1 ITIME(300),IOB(300),NBIAS(300) DOUBLE PRECISION DAL(IDIML) DOUBLE PRECISION DSAPR(300),DSAPO(300),DTLAG(300) CHARACTER CINSTR(300)*10 COMMON /BLOCKR/ IRECA,IRECE,IDATA,ITIMA,IDATE,ITIME,IOB,NBIAS, 1 DSAPR,DSAPO,DTLAG COMMON /BLOCKC/ CINSTR INTEGER*4 IDATL(20) DOUBLE PRECISION DTF(23),DATL(20),DGI(6) CHARACTER CUNIT(11)*8 COMMON /UNITS/ CUNIT,IC2 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DEFINE COEFFICIENTS OF 23 HOUR NUMERICAL DATA TEST FILTER, C C CF. WENZEL 1976, PAGE XVII. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DATA DTF/ -0.033261D0, 0.070770D0, 0.012603D0,-0.061752D0, 1-0.061916D0, 0.013218D0, 0.093693D0, 0.095321D0,-0.019683D0, 2-0.214064D0,-0.394930D0, 1.000000D0,-0.394930D0,-0.214064D0, 3-0.019683D0, 0.095321D0, 0.093693D0, 0.013218D0,-0.061916D0, 4-0.061752D0, 0.012603D0, 0.070770D0,-0.033261D0/ DATA MAXDE/20/,NFI/23/,NFI2/11/ CALL DATUM(IDATA(JB),ITIMA(JB),ITY,ITM,ITD,ITH) ITH=ITH+NFI2 IF(ITH.GE.24) CALL ETGREI(ITY,ITM,ITD,ITH) MM=NO-NFI+1 DC2=0.D0 NDLB=0 DO 290 J=1,MM,6 DO 270 N=1,6 ISTART=J+N-1 IEND=J+N+NFI-2 IF(IEND.GT.NO) GOTO 280 C CALL ETNFIL(DTF,NFI,NFI,NFI2,DAL,IDIML,ISTART,DFLOW,DFHIG) C IF(DABS(DFLOW).LE.DATLIM) GOTO 260 IF(NDLB.GE.MAXDE) GOTO 260 NDLB=NDLB+1 IDATL(NDLB)=J+N-1 DATL(NDLB)=DFLOW 260 DC2=DC2+DFLOW**2 NL=N 270 DGI(N)=DFLOW 280 IF(MOD(J,720).EQ.1) WRITE(IUN6,7036) JB WRITE(IUN6,7033) CINSTR(JB),ITY,ITM,ITD,ITH,(DGI(N),N=1,NL) ITH=ITH+6 IF(ITH.GE.24) CALL ETGREI(ITY,ITM,ITD,ITH) 290 CONTINUE DC2=DSQRT(DC2/DBLE(MM)) WRITE(IUN6,7037) JB,DC2,CUNIT(IC2) WRITE(*,7037) JB,DC2,CUNIT(IC2) 300 CONTINUE IF(NDLB.EQ.0) RETURN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C LIST OF DATA ERRORS. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC WRITE(IUN6,7044) NDLB,DATLIM,CUNIT(IC2),JB WRITE(*,7044) NDLB,DATLIM,CUNIT(IC2),JB DO 310 N=1,NDLB CALL DATUM(IDATA(JB),ITIMA(JB),ITY,ITM,ITD,ITH) ITH=ITH+IDATL(N)+10 IF(ITH.GE.24) CALL ETGREI(ITY,ITM,ITD,ITH) WRITE(IUN6,7051) ITY,ITM,ITD,ITH,DATL(N) 310 WRITE(*,7051) ITY,ITM,ITD,ITH,DATL(N) RETURN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Format statements: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 7033 FORMAT(A4,I5,3(1H.,I2),6F8.2) 7036 FORMAT(' Routine ETSDER version 930629 FORTRAN 77.'// 1 23H SEARCH FOR DATA ERRORS,5X,' BLOCK NO.',I4// 2 11H DATE IN UT,12X,2H+0,6X,2H+1,6X,2H+2,6X,2H+3,6X,2H+4,6X,2H+5/) 7037 FORMAT(//' OBSERVATION BLOCK NO. :',I10/ 1' RMS OF DATA ERRORS :',F10.3,2X,A8/) 7044 FORMAT(' *****',I5,' DATA ERRORS EXCEED GIVEN LIMIT OF',F10.3, 1 2X,A8,' IN BLOCK NO.',I4// 2' DATES OF DATA ERRORS:'/) 7051 FORMAT(1X,I4,1H.,I2,1H.,I2,1H.,I2,F10.3) END C SUBROUTINE GEOEXT(IUN6,DEXTIM) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Routine GEOEXT, version 930418 FORTRAN 77. C C C C === MS-DOS version for LAHEY-compiler =================== C C C C The routine GEOEXT computes the actual job time and writes C C the actual execution time on printer output unit IUN6. C C For the first call of routine GEOEXT, the actual jobtime will C C be computed (in secs since midnight) and stored. For the next C C call(s) of routine GEOEXT, the actual jobtime will be computed C C and the execution time (actual jobtime minus jobtime of the C C first call of routine GEOEXT) will be printed. C C C C Input parameter description: C C ---------------------------- C C C C IUN6... formatted printer unit. C C C C Output parameter description: C C ----------------------------- C C C C DEXTIM... actual jobtime in seconds (time elapsed from the C C first call of routine GEOEXT to the actual call of C C routine GEOEXT.), double precision. C C C C used routines: C C -------------- C C C C Subroutine GETTIM for WATFOR87 or IBM PROFESSIONAL FORTRAN C C compiler. C C C C Program creation: 790830 by H.-G. Wenzel, C C Geodaetisches Institut, C C Universitaet Karlsruhe, C C Englerstr. 7, C C D-76128 KARLSRUHE 1, C C Germany. C C Tel.: 0721-6082301. C C FAX: 0721-694552. C C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de C C Last Modification: 930418 by H.-G.Wenzel. C C**********************************************************************C IMPLICIT DOUBLE PRECISION (D) C MSFOR: INTEGER*2 IH,IM,IS,IS100 SAVE IZ1,IZ2,DTIME1 DATA IZ1/0/,IZ2/0/ IF(IZ1-IZ2) 6003,6001,6003 6001 CONTINUE IZ1=10000 C MSFOR: CALL GETTIM(IH,IM,IS,IS100) C MSFOR: DTIME1=DBLE(IS+IM*60+IH*3600)+0.01*FLOAT(IS100) C LAHEY: CALL TIMER(ITS100) DTIME1=DBLE(ITS100)/100.D0 C UNIX: DTIME1=DBLE(SECNDS(RDUMMY)) WRITE(IUN6,7001) WRITE(*,7001) DEXTIM=0.D0 RETURN 6003 CONTINUE C MSFOR: CALL GETTIM(IH,IM,IS,IS100) C MSFOR: DTIME2=DBLE(IS+IM*60+IH*3600)+0.01*FLOAT(IS100) C LAHEY: CALL TIMER(ITS100) DTIME2=DBLE(ITS100)/100.D0 C UNIX: DTIME2=DBLE(SECNDS(RDUMMY)) DEXTIM=DTIME2-DTIME1 WRITE(IUN6,7002) DEXTIM WRITE(*,7002) DEXTIM CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Format statements: C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 7001 FORMAT(' First call of routine GEOEXT, version 930418 FTN77.') 7002 FORMAT(/' Routine GEOEXT. Execution time=',F10.3,' sec'/) RETURN END
Become a Patron
Sponsor on GitHub
Donate via PayPal
Source on GitHub
Mailing list
Installed libraries
Wiki
Report an issue
How it works
Contact the author
CE on Mastodon
About the author
Statistics
Changelog
Version tree