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)
flang-trunk (flang-to-external-fc)
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 PREDICT C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Modified 2004.02.18 for Compaq Visual Fortran (B.Ducarme) ! C Program PREDICT, version 3.31 1997.03.03 Fortran 90. ! C ! C This file has 3707 records. ! C ! C === Version for MS-DOS using LAHEY LF90 Fortran compiler === ! C ! C === To run this program under UNIX, you have to modify === ! C routine GEOEXT. ! C ! C The program PREDICT computes model tides using different tidal ! C potential catalogues: ! C ! C Doodson (1921) with 378 waves, ! C Cartwright-Tayler-Edden (1973) with 505 waves, ! C Buellesfeld (1985) with 656 waves, ! C Tamura (1987) with 1200 waves, ! C Xi (1989) with 2933 waves, ! C Roosbeek (1996) with 6499 waves, ! C Hartmann and Wenzel (1995) with 12935 waves ! C ! C for a number of different tidal components using observed or ! C estimated (e.g. by a body tide and ocean tide model) tidal ! C parameters. ! C ! C Reference: ! C ---------- ! C ! C Wenzel, H.-G. (1996): The nanogal software: Earth tide data ! C processing package ETERNA 3.3. Bulletin d'Informations ! C Marees Terrestres vol. 124, 9425-9439, Bruxelles 1996. ! C ! C Disc file description: ! C ---------------------- ! C ! C project: Formatted file, on which the project name 'CPROJ' ! C has to be stored in the first record (8 characters ! C at maximum). The project name 'CPROJ' will be used ! C to define the print file and the output file. ! C This file has to be stored in the directory, from ! C which program PREDICT is executed. ! C 'CPROJ'.ini: Formatted unit, on which the input parameters ! C have to be stored before the execution of program ! C PREDICT. This file has to be stored in the ! C directory, from which program PREDICT is executed. ! C 'CPROJ'.prn: Formatted print file. This file will be written in ! C the directory, from which program PREDICT is ! C executed. ! C 'CPROJ'.prd: Formatted output file in ETERNA format. This file ! C will be written in the directory, from which ! C program PREDICT is executed. ! C doodsehw.dat: Formatted file, on which the Doodson (1921) tidal ! C potential catalogue has to be stored before the ! C execution of program PREDICT. ! C The path for this file is ! C /home/hwz/eterna34/commdat/doodsehw.dat. ! C doodseh2.uft: Unformatted file, on which the Doodson (1921) ! C tidal potential catalogue will be stored by the ! C first execution of program predict, if it does not ! C yet exist. The path for this file is ! C /home/hwz/eterna34/commdat/doodsehw.uft. ! C cted73hw.dat: Formatted file, on which the Cartwright and Tayler ! C (1971) and Cartwright and Edden (1973) tidal ! C potential catalogue has to be stored before the ! C execution of program PREDICT. ! C The path for this file is ! C /home/hwz/eterna34/commdat/cted73hw.dat. ! C cted73h2.uft: Unformatted file, on which the Cartwright and ! C Tayler (1971) and Cartwright and Edden (1973) ! C tidal potential catalogue will be stored by the ! C first execution of program predict, if it does not ! C yet exist. The path for this file is ! C /home/hwz/eterna34/commdat/cted73hw.uft. ! C buellehw.dat: Formatted file, on which the Buellesfeld (1985) ! C tidal potential catalogue has to be stored before ! C the execution of program PREDICT. ! C The path for this file is ! C /home/hwz/eterna34/commdat/buellehw.dat. ! C buelleh2.uft: Unformatted file, on which the Buellesfeld (1985) ! C tidal potential catalogue will be stored by the ! C first execution of program predict, if it does not ! C yet exist. The path for this file is ! C /home/hwz/eterna34/commdat/buellehw.uft. ! C tamurahw.dat: Formatted file, on which the Tamura (1987) ! C tidal potential catalogue has to be stored before ! C the execution of program PREDICT. ! C The path for this file is ! C /home/hwz/eterna34/commdat/tamurahw.dat. ! C tamurah2.uft: Unformatted file, on which the Tamura (1987) ! C tidal potential catalogue will be stored by the ! C first execution of program predict, if it does not ! C yet exist. The path for this file is ! C /home/hwz/eterna34/commdat/tamurah2.uft. ! C xi1989hw.dat: Formatted file, on which the Xi (1989) tidal ! C potential catalogue has to be stored before the ! C execution of program PREDICT. ! C The path for this file is ! C /home/hwz/eterna34/commdat/xi1989hw.dat. ! C xi1989h2.uft: Unformatted file, on which the Xi (1989) tidal ! C potential catalogue will be stored by the first ! C execution of program predict, if it does not yet ! C exist. The path for this file is ! C /home/hwz/eterna34/commdat/xi1989hw.uft. ! C ratgp95.dat: Formatted file, on which the Roosbeek (1986) tidal ! C potential catalogue has to be stored before the ! C execution of program PREDICT. ! C The path for this file is ! C /home/hwz/eterna34/commdat/ratgp95.dat. ! C ratgp952.uft: Unformatted file, on which the Roosbeek (1986) ! C tidal potential catalogue will be stored by the ! C first execution of program predict, if it does not ! C yet exist. The path for this file is ! C /home/hwz/eterna34/commdat/ratgp95.uft. ! C hw95s.dat: Formatted file, on which the Hartmann and Wenzel ! C (1995) tidal potential catalogue has to be stored ! C before the execution of program PREDICT. ! C The path for this file is ! C /home/hwz/eterna34/commdat/hw95.dat. ! C hw95s2.uft: Unformatted file, on which the Hartmann and Wenzel ! C (1995) tidal potential catalogue will be stored by ! C the first execution of program PREDICT, if it does ! C not yet exist. The path for this file is ! C /home/hwz/eterna34/commdat/hw95.uft. ! C etpolut1.dat: Formatted file, on which the pole coordinates and ! C DUT1 corrections have to be stored before the ! C execution of program PREDICT. ! C The path for this file is ! C /home/hwz/eterna34/commdat/etpolut1.dat. ! C etpolut2.uft: Unformatted direct access file, on which the pole ! C coordinates and DUT1 corrections will be stored by ! C the first execution of program PREDICT, if it does ! C not yet exist. The path for this file is ! C /home/hwz/eterna34/commdat/etpolut2.uft. ! C ! C Used routines: ! C -------------- ! C ! C ETASTN: computes astronomical elements. ! C PREDIN: reads control parameters. ! C ETDDTA: reads tabel of DDT = TDT - UTC. ! C ETDDTB: interpolates DDT = TDT - UTC from table. ! C ETGCON: computes geodetic coefficients. ! C ETGREN: computes date from Julian date ! C ETJULN: computes JULIAN date. ! C ETLEGN: computes fully normalized Legendre spherical harmonics. ! C ETLOVE: computes elastic parameters from Wahr-Dehant model. ! C ETPOLC: computes DUT1 ! C ETPHAS: computes the phases and frequencies of the tidal waves. ! C ETPOTS: computes amplitudes, frequencies and phases of tidal ! C waves. ! C GEOEXT: computes JOBTIME. ! C ! C Numerical accuracy: ! C ------------------- ! C ! C The program has been tested on IBM-At compatible computers under ! C MS-DOS operating system with different compilers using DOUBLE ! C PRECISION for all variables (15 digits) and on a SUN SPARC2 ! C under UNIX operating system with SUN F77 compiler, and gave ! C identical results to 0.0001 nm/s**2. ! C ! C Execution time: ! C --------------- ! C ! C The CPU execution time depends mainly on the number of waves of ! C the choosen tidal potential development, and the number of model ! C tide values to be computed. ! C The execution on a PENTIUM 100 MHz using LF90 compiler for 8760 ! C hourly samples including pole tide and LOD tide (parameter file ! C KAHW9501.INI) are ! C ! C catalogue threshold nwave rms error ex. time ! C [nm/s**2] [sec] ! C ! C Tamura (1987) 2.D-06 1200 0.070 11.86 ! C ! C Hartmann + Wenzel (1995) 1.D-01 9 88.40 3.90 ! C Hartmann + Wenzel (1995) 1.D-02 42 14.40 4.06 ! C Hartmann + Wenzel (1995) 1.D-03 155 2.25 4.94 ! C Hartmann + Wenzel (1995) 1.D-04 434 0.44 7.20 ! C Hartmann + Wenzel (1995) 1.D-05 1248 0.068 13.51 ! C Hartmann + Wenzel (1995) 1.D-06 3268 0.011 33.01 ! C Hartmann + Wenzel (1995) 1.D-07 7761 0.002 83.82 ! C Hartmann + Wenzel (1995) 1.D-08 11462 0.001 132.86 ! C Hartmann + Wenzel (1995) 1.D-09 12000 0.001 139.95 ! C Hartmann + Wenzel (1995) 1.D-10 12011 0.001 140.11 ! C ! C Program creation: 1973.06.23 by Hans-Georg Wenzel, ! C Black Forest Observatory, ! C Universitaet Karlsruhe, ! C Englerstr. 7, ! C D-76128 KARLSRUHE, ! C Germany. ! C Tel.: 0721-6082307, ! C FAX: 0721-694552. ! C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de ! C Last modification: 1997.03.03 by Hans-Georg Wenzel. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IMPLICIT DOUBLE PRECISION (D) CHARACTER CUNIT(11)*8,COMPON(11)*24,CVERS*11,C88*8,C99*9 CHARACTER CHEAD(10)*64 CHARACTER CPROJ*8,CFINI*13,CFPRN*13,CFOUT*13 DIMENSION DGI(6) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C The following DIMENSION statements are concerning the number of ! C output channels: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! PARAMETER (MAXCHAN=4) CHARACTER CHANNEL(MAXCHAN)*10 DIMENSION DCOUT(MAXCHAN),DZERO(MAXCHAN) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C The following DIMENSION statement is concerning the tabel of ! C differences DDT = TDT - UTC: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DOUBLE PRECISION DDTTAB(3,300) COMMON /DDT/ DDTTAB,NDDTAB C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C The following dimension statements are concerning the number of ! C wavegroups to be used, which is restricted to 85 in the current ! C program version (parameter MAXWG). ! C ! C The following list of wavegroups may be used for the analyis and ! C prediction of tides (frequencies in cycle per day referring to ! C epoch J2000): ! C ! C no. from to frequency of name of ! C freqency frequency main wave main wave ! C ! C 1 0.000000 0.001369 0.000000 long ! C 2 0.000133 0.004107 0.002738 SA ! C 3 0.004108 0.020884 0.005476 SSA ! C 4 0.020885 0.054747 0.036292 MM ! C 5 0.054748 0.091348 0.073202 MF ! C 6 0.091349 0.501369 0.109494 MTM ! C 7 0.501370 0.911390 0.893244 Q1 ! C 8 0.911391 0.947991 0.929536 O1 ! C 9 0.947992 0.981854 0.966446 M1 ! C 10 0.981855 0.998631 0.997262 P1 ! C 11 0.998632 1.001369 1.000000 S1 ! C 12 1.001370 1.004107 1.002738 K1 ! C 13 1.004108 1.006845 1.005476 PSI1 ! C 14 1.006846 1.023622 1.008214 PHI1 ! C 15 1.023623 1.057485 1.039030 J1 ! C 16 1.057486 1.470243 1.075940 OO1 ! C 17 1.470244 1.880264 1.864547 2N2 ! C 18 1.880265 1.914128 1.895982 N2 ! C 19 1.914129 1.950419 1.932274 M2 ! C 20 1.950420 1.984282 1.968565 L2 ! C 21 1.984283 2.002736 2.000000 S2 ! C 22 2.002737 2.451943 2.005476 K2 ! C 23 2.451944 7.000000 2.898410 M3M6 ! C ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! PARAMETER (MAXWG=85) DIMENSION DFRA(MAXWG),DFRE(MAXWG),NA(MAXWG),NE(MAXWG), 1 DG0(MAXWG),DPHI0(MAXWG),DAM(MAXWG),DBOD(MAXWG) CHARACTER CNSY(MAXWG)*4 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C The following common blocks are used to transfer the control ! C parameters from routine PREDIN to main program. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! COMMON /CONTROL3/ DLAT,DLON,DH,DGRAV,DAZ,DFRA,DFRE,DG0,DPHI0, 1 DATLIM,DAMIN,DPOLTC,DLODTC,IDTSEC CHARACTER CINST*10 PARAMETER (MAXNF=8) INTEGER IREG(MAXNF) CHARACTER CFY1(MAXNF)*10,CFY2(MAXNF)*10 COMMON /CONTROL4/ IC,IR,ITY,ITM,ITD,ITH,IDA,KFILT,IPROBS, 1 IPRLF,IMODEL,IRIGID,IHANN,IQUICK,ISPANH,NGR,NF, 2 IREG,CFY1,CFY2,CINST,CNSY,CHEAD C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C The following dimension statements are concerning the number of ! C waves of the tidal potential catalogue, which is 13 000 in the ! C current program version (parameter MAXNW). ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! PARAMETER (MAXNW=12935) DOUBLE PRECISION DC0(MAXNW),DS0(MAXNW),DDC(MAXNW),DDS(MAXNW) COMMON /TIDWAVE/ NW,IWNR(12935),IAARG(12935,12),DX0(12935), 1 DX1(12935),DY0(12935),DY1(12935),DTHPH(12935),DTHFR(12935), 2 DBODY(12935) COMMON /UNITS/ CUNIT,IC2 COMMON /CONST/ DPI,DPI2,DRAD,DRO DATA IUN14/14/,IUN15/15/,IUN16/16/,IUN23/23/,IUN24/24/,IUN27/27/ DATA IUN30/30/,IUN31/31/ DATA DZERO/4*0.D0/ DATA CVERS/'3.31 970303'/,C88/'88888888'/,C99/'999999999'/ DATA COMPON/'Potential ' ,'Gravity ', 1'Tilt ','Vertical displacement ', 3'Horizontal displacement ','Vertical strain ', 5'Horizontal strain ','Aereal strain ', 7'Shear strain ','Volume strain ', 9'Ocean tide '/ C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Read project name: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! OPEN(UNIT=IUN15,FILE='project') READ(IUN15,17001) CPROJ CLOSE(IUN15) CFINI=CPROJ//'.ini' CFPRN=CPROJ//'.prn' CFOUT=CPROJ//'.prd' C OPEN(IUN15, FILE=CFINI,FORM='FORMATTED') OPEN(IUN16, FILE=CFPRN,FORM='FORMATTED') OPEN(IUN23, FILE=CFOUT,FORM='FORMATTED') WRITE(IUN16,17002) CVERS,CPROJ WRITE(*,17002) CVERS,CPROJ IRESET=1 CALL GEOEXT(IUN16,IRESET,DEXTIM,DEXTOT) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Store array of differences DDT = TDT - UTC: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IPRINT=0 CALL ETDDTA(IUN16,IUN27,IPRINT) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Read control parameters: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IPRINT=1 CALL PREDIN(IUN15,IUN16,IPRINT) DDTH=DBLE(IDTSEC)/3600.D0 DDTD=DDTH/24.D0 DTH=0.D0 IF(IRIGID.EQ.1) WRITE(IUN16,17009) WRITE(IUN23,17018) CVERS C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Define output channels: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IPOLTC=0 ILODTC=0 NC=1 CHANNEL(1)='pred.sig. ' IF(DABS(DPOLTC).GT.1.D-6) IPOLTC=1 IF(DABS(DLODTC).GT.1.D-6) ILODTC=1 IF(IPOLTC.EQ.1.OR.ILODTC.EQ.1) THEN NC=NC+1 CHANNEL(NC)='pred.tide ' ENDIF IF(IPOLTC.EQ.1) THEN NC=NC+1 CHANNEL(NC)='pole tide ' ENDIF IF(ILODTC.EQ.1) THEN NC=NC+1 CHANNEL(NC)='lod tide ' ENDIF IPRINT=1 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Initialize direct access file for polecoordinates and DUT1: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C HP-UNIX: OPEN(UNIT=IUN30,FILE='../commdat/etpolut1.dat', C HP-UNIX: 1 FORM='FORMATTED',STATUS='OLD') C MS-DOS: OPEN(UNIT=IUN30,FILE='/home/hwz/eterna34/commdat/etpolut1.dat', 1 FORM='FORMATTED',STATUS='OLD') DCLAT=DCOS(DLAT*DRAD) DSLAT=DSIN(DLAT*DRAD) DCLON=DCOS(DLON*DRAD) DSLON=DSIN(DLON*DRAD) IPRINT=1 CALL ETJULN(IUN16,ITY,ITM,ITD,DTH,DJULD) CALL ETPOLC(IUN16,IUN30,IUN31,IPRINT,DJULD,DCLAT,DSLAT, 1 DCLON,DSLON,DPOLX,DPOLY,DUT1,DTAI,DLOD,DGPOL,DGPOLP,DGLOD,NERR) CLOSE(IUN30) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Compute amplitudes, frequencies and phases of the tidal waves. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IPRINT=1 CALL ETPOTS(IUN14,IUN16,IUN24,IPRINT,IMODEL,DLAT,DLON,DH,DGRAV, 1 DAZ,IC,DJULD,DAMIN) IC2=IC+2 ITH=DTH C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Print output channel table: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DO 100 JC=1,NC WRITE(IUN16,17003) JC,CHANNEL(JC),CUNIT(IC2) 100 WRITE(IUN23,17003) JC,CHANNEL(JC),CUNIT(IC2) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Print alphanumeric comment: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! WRITE(IUN16,17010) CVERS DO 900 I=1,10 900 WRITE(IUN16,17012) CHEAD(I) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Check wave groups and observed tidal parameters: ! C DG0: amplitude factor. ! C DPHI0... phase lead in degree. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! WRITE(IUN16,17013) CUNIT(IC2) WRITE(IUN23,17013) CUNIT(IC2) DO 910 IG=1,NGR C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Convert frequencies from cpd to rad per hour: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DFRA(IG)=DFRA(IG)*15.D0*DRAD DFRE(IG)=DFRE(IG)*15.D0*DRAD DO 930 IW=1,NW IF(DTHFR(IW).LT.DFRA(IG)-1.D-10) NA(IG)=IW+1 IF(DTHFR(IW).LT.DFRE(IG)+1.D-10) NE(IG)=IW 930 CONTINUE IF(NA(IG).EQ.0) NA(IG)=1 NAK=NA(IG) NEK=NE(IG) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Search for main wave of the group: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DAM(IG)=0.D0 DO 970 IW=NAK,NEK IF(IRIGID.EQ.1) DBODY(IW)=1.D0 DTHAM=DSQRT(DX0(IW)**2+DY0(IW)**2) IF(DTHAM.LT.DAM(IG)) GOTO 970 DAM(IG)=DTHAM DBOD(IG)=DBODY(IW) 970 CONTINUE 910 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Check last group: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 980 IF(NE(NGR).GE.NA(NGR)) GOTO 990 NGR=NGR-1 GOTO 980 990 CONTINUE DO 995 IG=1,NGR NANR=IWNR(NA(IG)) NENR=IWNR(NE(IG)) WRITE(IUN16,17015) IG,NANR,NENR,DG0(IG),DPHI0(IG),CNSY(IG), 1 DAM(IG),DBOD(IG) WRITE(IUN23,17015) IG,NANR,NENR,DG0(IG),DPHI0(IG),CNSY(IG), 1 DAM(IG),DBOD(IG) DPHI0(IG)=DPHI0(IG)*DRAD 995 CONTINUE DO 996 IG=1,NGR DFAC=DG0(IG)/DBOD(IG) DO 996 IW=NA(IG),NE(IG) IF(IRIGID.EQ.1) DBODY(IW)=1.D0 DX0(IW)=DX0(IW)*DFAC*DBODY(IW) DX1(IW)=DX1(IW)*DFAC*DBODY(IW) DY0(IW)=DY0(IW)*DFAC*DBODY(IW) DY1(IW)=DY1(IW)*DFAC*DBODY(IW) 996 CONTINUE CALL GEOEXT(IUN16,IRESET,DEXTIM,DEXTOT) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Print hourly model tides with format 6F13.6: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1000 WRITE(IUN16,17016) CVERS,COMPON(IC2),CUNIT(IC2) WRITE(IUN23,17019) WRITE(IUN23,17020) (DZERO(J),J=1,NC) ITMIN=0 ITSEC=0 DDAT=DBLE(ISPANH)/DDTH NDAT=DDAT CALL ETJULN(IUN16,ITY,ITM,ITD,DTH,DJULD0) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Loop over NDAT samples: ! C DTLIM is the time interval in hours for updatinmg the phases. ! C Depending on amplitude threshold DAMIN, the phases will be ! C updated ! C at each midnight for DAMIN <= 1.D-8 m**2/s**2 ! C at monthly interval for 1.D-8 < DAMIN <= 1.D-6 m**2/s**2) ! C at yearly interval for 1.D-6 < DAMIN. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DT=1.D99 DTLIM=1.D0 IF(DAMIN.GT.1.D-8) DTLIM=720.D0 IF(DAMIN.GT.1.D-6) DTLIM=8760.D0 DO 1090 I=1,NDAT,6 DO 1020 J=1,6 DJULD=DJULD0+DBLE(I+J-2)*DDTD DT2000=(DJULD-2451544.D0)/36525.0D0 CALL ETGREN(IUN16,DJULD,ITY,ITM,ITD,DTH,NERR) IF(DT.LT.DTLIM) GOTO 1033 IF(DTLIM.LT.10.D0.AND.DTH.GT.0.0001D0) GOTO 1033 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Recompute phases at interval DTLIM: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IPRINT=0 CALL ETPHAS(IUN16,IPRINT,IMODEL,DLON,DJULD) DO 1025 IG=1,NGR DO 1025 IW=NA(IG),NE(IG) DTHPH(IW)=DTHPH(IW)+DPHI0(IG) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Prepare arrays for recursion algorithm: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DC0(IW)=DCOS(DTHPH(IW)) DS0(IW)=DSIN(DTHPH(IW)) DDC(IW)=DCOS(DTHFR(IW)*DDTH) 1025 DDS(IW)=DSIN(DTHFR(IW)*DDTH) DT=0.D0 1033 CONTINUE DGT=0.D0 DO 1030 IG=1,NGR DO 1030 IW=NA(IG),NE(IG) DGT=DGT+(DX0(IW)+DT2000*DX1(IW))*DC0(IW)+ 1 (DY0(IW)+DT2000*DY1(IW))*DS0(IW) DUMMY =DC0(IW)*DDC(IW)-DS0(IW)*DDS(IW) DS0(IW)=DS0(IW)*DDC(IW)+DC0(IW)*DDS(IW) 1030 DC0(IW)=DUMMY C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Compute pole tide correction and length of day tide correction: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DPOLT=0.D0 DLODT=0.D0 IF(IC.NE.0) GOTO 1040 IF(IPOLTC.EQ.1.OR.ILODTC.EQ.1) THEN CALL ETPOLC(IUN16,IUN30,IUN31,IPRINT,DJULD,DCLAT,DSLAT,DCLON, 1 DSLON,DPOLX,DPOLY,DUT1,DTAI,DLOD,DGPOL,DGPOLP,DGLOD,IKENN) DPOLT=DGPOL*DPOLTC DLODT=DGLOD*DLODTC ENDIF 1040 CONTINUE DCOUT(1)=DGT+DPOLT+DLODT DCOUT(2)=DGT JC=2 IF(IPOLTC.EQ.1) THEN JC=JC+1 DCOUT(JC)=DPOLT ENDIF IF(ILODTC.EQ.1) THEN JC=JC+1 DCOUT(JC)=DLODT ENDIF DGI(J)=DCOUT(1) IDAT=ITY*10000+ITM*100+ITD ITH=DTH DTMIN=(DTH-DBLE(ITH))*60.D0 ITMIN=INT(DTMIN+0.5D0) IF(ITMIN.GE.60) THEN ITMIN=0 ITH=ITH+1 ENDIF DTSEC=DTH*3600.D0-DBLE(ITH)*3600.D0-DBLE(ITMIN)*60.D0 ITSEC=INT(DTSEC+0.5D0) ITIM=ITH*10000+ITMIN*100+ITSEC WRITE(IUN23,17021) IDAT,ITIM,(DCOUT(JC),JC=1,NC) IF(ITIM.EQ.0) WRITE(*,17022) IDAT,ITIM,(DCOUT(JC),JC=1,NC) DT=DT+DDTH 1020 CONTINUE DJULD=DJULD0+DBLE(I-1)*DDTD CALL ETGREN(IUN16,DJULD,ITY,ITM,ITD,DTH,NERR) IDAT=ITY*10000+ITM*100+ITD ITH=DTH DTMIN=(DTH-DBLE(ITH))*60.D0 ITMIN=INT(DTMIN+0.5D0) IF(ITMIN.GE.60) THEN ITMIN=0 ITH=ITH+1 ENDIF DTSEC=DTH*3600.D0-DBLE(ITH)*3600.D0-DBLE(ITMIN)*60.D0 ITSEC=INT(DTSEC+0.5D0) ITIM=ITH*10000+ITMIN*100+ITSEC WRITE(IUN16,17017) IDAT,ITIM,(DGI(J),J=1,6) 1090 CONTINUE WRITE(IUN23,17001) C99 WRITE(IUN23,17001) C88 CALL GEOEXT(IUN16,IRESET,DEXTIM,DEXTOT) WRITE(IUN16,17030) CPROJ,CFPRN,CFOUT,DEXTIM WRITE(*,17030) CPROJ,CFPRN,CFOUT,DEXTIM C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Format statements: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 17001 FORMAT(A8) 17002 FORMAT( 1' ******************************************************'/ 2' * *'/ 3' * Program PREDICT, version ',A11,' Fortran 90. *'/ 4' * *'/ 5' * Prediction of earthtide signals. *'/ 6' * for project ',A8, ' *'/ 7' * *'/ 8' * The Black Forest Observatory Schiltach *'/ 9' * wishes you much success when using PREDICT. *'/ *' * *'/ 1' ******************************************************'/) 17003 FORMAT(6x,'output channel no. ',I5,2X,A10,2X,A10) 17009 FORMAT(//' ***** Parameter IRIGID = 1 has been input. '/ 1' ***** IRIGID =1 should only be used for comparison with model'/ 2' ***** tides computed from ephemeris programs. For real world '/ 3' ***** computations, use always IRIGID=0.'/) 17010 FORMAT(//6x,'Program PREDICT, version ',A11,'Fortran 90'//) 17012 FORMAT(1X,A64) 17013 FORMAT(/// 1 6x,'Wave groups and observed tidal parameters'// 2 6x,' no. from to ampl.fac. phase lead ampl. WD body '/ 3 6x,' [deg] [',A10,']'/) 17015 FORMAT(6x,3I6,2F10.4,1X,A8,2F10.4) 17016 FORMAT(//// 1 6x,'Program PREDICT, version ',A11,' Fortran 90'// 2 6x,'Component ',A24,' IN ',1X,A8//6X,'Date in UT'//) 17017 FORMAT(5X,I8,1X,I6,6F10.3) 17018 FORMAT(' Model tides from program PREDICT, version',A11/) 17019 FORMAT(/'C******************************************************'/ 1'PREDICT 1.0000 1.0000 0.000 3', 2'PREDICT') 17020 FORMAT('77777777',7X,4F10.3) 17021 FORMAT(I8,1X,I6,6F10.3) 17022 FORMAT(I9,1X,I2,6F10.3) 17030 FORMAT(/ 1' **********************************************'/ 2' * Program PREDICT finished the execution *'/ 3' * for project ',A8, ' *'/ 4' * (Hopefully it was successfull). *'/ 5' * Print file is: ',A13,' *'/ 6' * Output file is: ',A13,' *'/ 7' **********************************************'// 6' Execution time: ',F10.3,' seconds'/) END C BLOCK DATA C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C ! C BLOCK DATA for program PREDICT, version 1996.05.25 Fortran 90. ! C ! C Routine creation: 1993.03.29 by Hans-Georg Wenzel, ! C Black Forest Observatory, ! C Universitaet Karlsruhe, ! C Englerstr. 7, ! C D-76128 KARLSRUHE 1, ! C Germany. ! C Tel: 0049-721-6082307, ! C FAX: 0049-721-694552. ! C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de ! C Last modification: 1996.05.25 by Hans-Georg Wenzel. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IMPLICIT DOUBLE PRECISION (D) CHARACTER CUNIT(11)*8 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C COMMON /CONST/: ! C DPI... 3.1415.... DPI2... 2.D0*DPI ! C DRAD... DPI/180.D0 DRO... 180.D0/DPI ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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 ETASTN(IUN16,IPRINT,IMODEL,DLON,DJULD,DUT1,DAS,DASP, 1 DDT) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C ! C Routine ETASTN, version 1996.05.25 Fortran 90. ! C ! C The routine ETASTN computes the astronomical elements for ! C different tidal potential catalogues at a specific epoch, given ! C in UTC. The formulas for the astronomical elements have been ! C taken from Tamura (1987) and Simon et al. (1994). ! C ! C Reference: ! C ---------- ! C ! C Simon, J.L., P. Bretagnon, J. Chapront, M. Chapront-Touze, ! C G. Francou and J. Laskar (1994): Numerical expressions for ! C precession formulae and mean elements for the Moon and the ! C planets. Astronomy and Atsrohysics, vo. 282, 663-683, 1994. ! C Tamura, Y. (1987): A harmonic development of the tide ! C generating potential. Bulletin d'Informations Marees ! C Terrestres vol. 99, 6813-68755, Bruxelles 1987. ! C ! C All variables with D as first character are double precision. ! C ! C Input parameter description: ! C ---------------------------- ! C ! C IUN16: Unit number of formatted printout file. ! C IPRINT: Printout parameter. For IPRINT=0, nothing will ! C be printed on unit IUN16. ! C IMODEL: Parameter describing the tidal potential catalogue. ! C IMODEL = 1: Doodson (1921) catalogue. ! C IMODEL = 2: Cartwright et al. (1973) catalogue. ! C IMODEL = 3: Buellesfeld (1985) catalogue. ! C IMODEL = 4: Tamura (1987) catalogue. ! C IMODEL = 5: Xi (1989) catalogue. ! C IMODEL = 6: Roosbeek (1996) catalogue. ! C IMODEL = 7: Hartmann and Wenzel (1995) catalogue. ! C For IMODEL = 1...5, arguments are computed from ! C Tamura (1987) formulas. For IMODEL = 6 and 7, ! C arguments are computed from Simon et al. (1994) ! C formulas. ! C DJULD: Julian date of the epoch in UTC. ! C ! C Output parameter description: ! C ----------------------------- ! C ! C DAS(1): Mean local Moontime in degree. ! C DAS(2): Mean longitude of the Moon in degree. ! C DAS(3): Mean longitude of the Sun in degree. ! C DAS(4): Mean longitude of the perigee of the Moon's orbit ! C in degree. ! C DAS(5): Negative mean longitude of the ascending node of ! C the Moon's orbit in degree. ! C DAS(6): Mean longitude of the perigee of the Suns's orbit ! C in degree. ! C DAS(7): Mean longitude of the Mercury in degree. ! C DAS(8): Mean longitude of the Venus in degree. ! C DAS(9): Mean longitude of the Mars in degree. ! C DAS(10): Mean longitude of the Jupiter in degree. ! C DAS(11): Mean longitude of the Saturn in degree. ! C ! C DASP(1...11): Time derivatives of the corresponding variables ! C DAS in degree per hour. ! C ! C Used routines: ! C -------------- ! C ETDDTB: interpolates DDT = DTD - UTC from table. ! C ! C Routine creation: 1994.07.30 by Hans-Georg Wenzel, ! C Geodaetisches Institut, ! C Universitaet Karlsruhe, ! C Englerstr. 7, ! C D-76128 KARLSRUHE 1, ! C Germany. ! C Tel.: 0721-6082307, ! C FAX: 0721-694552. ! C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de ! C Last modification: 1996.05.25 by Hans-Georg Wenzel. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IMPLICIT DOUBLE PRECISION (D) IMPLICIT INTEGER (I-N) DOUBLE PRECISION DAS(11),DASP(11) SAVE DATA DRAD/0.174532925197721D-001/ D1MD=1.D0/(365250.D0*24.D0) DMJD=DJULD-2400000.5D0 IMJD=DMJD DTH=(DMJD-DBLE(IMJD))*24.D0 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Compute Universal Time epoch DTUT in Julian Centuries referring ! C to J2000: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DTUT=(DJULD-2451545.0D0)/36525.D0 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Correct DTH to UT1: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DTH=DTH+DUT1/3600.D0 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Compute epoch DT in Julian Centuries TDB referring to J2000 ! C (1. January 2000 12 h.): ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DT=(DMJD-51544.5D0)/36525.0D0 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Correct time from UTC to TDT: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CALL ETDDTB(IUN16,IPRINT,DJULD,DDT) DT=DT+DDT/3155760000.D0 IF(IPRINT.GT.0) WRITE(IUN16,17001) DMJD DT2=DT*DT DTC1=DT DTC2=DTC1*DTC1 DTC3=DTC2*DTC1 DTC4=DTC3*DTC1 DTM1=DT/10.D0 DTM2=DTM1*DTM1 DTM3=DTM2*DTM1 DTM4=DTM3*DTM1 DTM5=DTM4*DTM1 DTM6=DTM5*DTM1 IF(IMODEL.GE.6) GOTO 2000 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Compute astronomical elements from TAMURA's 1987 formulas: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DTUT2=DTUT*DTUT DTUT3=DTUT2*DTUT DAL=280.4606184D0 + 36000.7700536D0*DTUT + 0.00038793D0*DTUT2 1 -0.0000000258D0*DTUT3 DALP=(36000.7700536D0 +2.0D0*0.00038793D0*DTUT 1 -3.0D0*0.0000000258D0*DTUT2)/(24.0D0*36525.D0) DS=218.316656D0+481267.881342D0*DT-0.001330D0*DT2 DSP=(481267.881342D0-2.0D0*0.001330D0*DT)/(24.D0*36525.0D0) DH=280.466449D0+36000.769822D0*DT+0.0003036D0*DT2 DHP=(36000.769822D0+2.0D0*0.0003036D0*DT)/(24.D0*36525.0D0) DDS=0.0040D0*DCOS((29.D0+133.0D0*DT)*DRAD) DDSP=(-0.0040D0*133.0D0*DRAD*DSIN((29.D0+133.0D0*DT)*DRAD))/ 1 (24.0D0*36525.0D0) DDH=0.0018D0*DCOS((159.D0+19.D0*DT)*DRAD) DDHP=(-0.0018D0*19.0D0*DRAD*DSIN((159.D0+19.D0*DT)*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*DT -0.010324D0*DT2 DAS(5)=234.955444D0 +1934.136185D0*DT -0.002076D0*DT2 DAS(6)=282.937348D0 + 1.719533D0*DT +0.0004597D0*DT2 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Compute the speeds in degree per hour: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DASP(1)=DALP-DSP+15.0D0 DASP(2)=DSP+DDSP DASP(3)=DHP+DDHP DASP(4)=(4069.013711D0-2.0D0*0.010324D0*DT)/(24.0D0*36525.0D0) DASP(5)=(1934.136185D0-2.0D0*0.002076D0*DT)/(24.0D0*36525.0D0) DASP(6)=(1.719533D0+2.0D0*0.0004597D0*DT)/(24.0D0*36525.0D0) GOTO 3000 2000 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Mean longitude of the Moon (from Simon et al. 1994): ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DS =218.3166456300D0+481267.8811957500D0*DTC1 2 -0.0014663889D0*DTC2 3 +0.0000018514D0*DTC3 4 -0.0000000153D0*DTC4 DSP=(+481267.8811957500D0 2 -2.D0*0.0014663889D0*DTC1 3 +3.D0*0.0000018514D0*DTC2 4 -4.D0*0.0000000153D0*DTC3)/(36525.D0*24.D0) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Mean longitude of the Sun (from Simon et al. 1994): ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DH=280.46645016D0+ 360007.6974880556D0*DTM1 2 +0.0303222222D0*DTM2 3 +0.0000200000D0*DTM3 4 -0.0000653611D0*DTM4 DHP= (360007.6974880556D0 2 +2.D0*0.0303222222D0*DTM1 3 +3.D0*0.0000200000D0*DTM2 4 -4.D0*0.0000653611D0*DTM3)*D1MD DAS(1) =DH -DS +DLON+DTH*15.0D0 DASP(1)=DHP-DSP+15.0D0 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Modification for Roosbeek (1996) tidal potential catalogue: ! C This modification has been programmed by Roosbeek himself. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(IMODEL.EQ.6) THEN DGMST=280.460618375D0+360007.700536D0*DTM1 2 +0.038793333333D0*DTM2 3 -0.000025833333D0*DTM3 DGMSTP= (360007.700536D0 2 +2.D0*0.038793333333D0*DTM1 3 -3.D0*0.000025833333D0*DTM2)*D1MD DAS(1) =DGMST-DS+DLON+DTH*15.D0 DASP(1)=DGMSTP-DSP+15.D0 ENDIF C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C This correction is necessary because for the determination of ! C the HW95 tidal potential catalogue the difference DDT=TDT-UTC ! C has been neglected. If the GMST would have been computed with ! C with the correct DDT, the effect in GMST would be 1.0027*DDT. ! C This effect is corrected below. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DAS(1)=DAS(1)-0.0027D0*DDT*15.D0/3600.D0 DAS(2) =DS DASP(2)=DSP DAS(3) =DH DASP(3)=DHP C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Mean longitude of lunar perigee (from Simon et al. 1994): ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DAS(4)= 83.35324312D0+40690.1363525000D0*DTM1 2 -1.0321722222D0*DTM2 3 -0.0124916667D0*DTM3 4 +0.0005263333D0*DTM4 DASP(4)= (+40690.1363525000D0 2 -2.D0*1.0321722222D0*DTM1 3 -3.D0*0.0124916667D0*DTM2 4 +4.D0*0.0005263333D0*DTM3)*D1MD C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Negative mean longitude of the ascending node of the Moon ! C in degree (from Simon et al. 1994): ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DAS(5)=234.95544499D0+19341.3626197222D0*DTM1 2 -0.2075611111D0*DTM2 3 -0.0021394444D0*DTM3 4 +0.0001649722D0*DTM4 DASP(5)= (+19341.3626197222D0 2 -2.D0*0.2075611111D0*DTM1 3 -3.D0*0.0021394444D0*DTM2 4 +4.D0*0.0001649722D0*DTM3)*D1MD C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Mean longitude of solar perigee computed from ! C argument no. 2 - D -l': ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DAS(6)=282.93734098D0 +17.1945766666D0*DTM1 1 +0.0456888889D0*DTM2 2 -0.0000177778D0*DTM3 2 -0.0000334444D0*DTM4 DASP(6)= (+17.1945766666D0 1 +2.D0*0.0456888889D0*DTM1 2 -3.D0*0.0000177778D0*DTM2 2 -4.D0*0.0000334444D0*DTM3)*D1MD 3000 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Longitudes of the planets from Simon et al. 1994: ! C Mercury = 7, Venus = 8, Mars = 9, Jupiter = 10, Saturn = 11. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DAS( 7)=252.25090552D0+1494740.7217223248D0*DTM1 2 +0.0303498417D0*DTM2 3 +0.0000181167D0*DTM3 4 -0.0000652778D0*DTM4 5 -0.0000004972D0*DTM5 6 +0.0000000556D0*DTM6 DASP( 7)= (+1494740.7217223248D0 2 +2.D0*0.0303498417D0*DTM1 3 +3.D0*0.0000181167D0*DTM2 4 -4.D0*0.0000652778D0*DTM3 5 -5.D0*0.0000004972D0*DTM4 6 +6.D0*0.0000000556D0*DTM5)*D1MD DAS( 8)=181.97980085D0+ 585192.1295333027D0*DTM1 2 +0.0310139472D0*DTM2 3 +0.0000149111D0*DTM3 4 -0.0000653222D0*DTM4 5 -0.0000004972D0*DTM5 6 +0.0000000556D0*DTM6 DASP( 8)= (+585192.1295333027D0 2 +2.D0*0.0310139472D0*DTM1 3 +3.D0*0.0000149111D0*DTM2 4 -4.D0*0.0000653222D0*DTM3 5 -5.D0*0.0000004972D0*DTM4 6 +6.D0*0.0000000556D0*DTM5)*D1MD DAS( 9)=355.43299958D0+ 191416.9637029695D0*DTM1 2 +0.0310518722D0*DTM2 3 +0.0000156222D0*DTM3 4 -0.0000653222D0*DTM4 5 -0.0000005000D0*DTM5 6 +0.0000000556D0*DTM6 DASP( 9)= (+191416.9637029695D0 2 +2.D0*0.0310518722D0*DTM1 3 +3.D0*0.0000156222D0*DTM2 4 -4.D0*0.0000653222D0*DTM3 5 -5.D0*0.0000005000D0*DTM4 6 +6.D0*0.0000000556D0*DTM5)*D1MD DAS(10)= 34.35151874D0+ 30363.0277484806D0*DTM1 2 +0.0223297222D0*DTM2 3 +0.0000370194D0*DTM3 4 -0.0000523611D0*DTM4 5 +0.0000011417D0*DTM5 6 -0.0000000389D0*DTM6 DASP(10)= (+30363.0277484806D0 2 +2.D0*0.0223297222D0*DTM1 3 +3.D0*0.0000370194D0*DTM2 4 -4.D0*0.0000523611D0*DTM3 5 +5.D0*0.0000011417D0*DTM4 6 -6.D0*0.0000000389D0*DTM5)*D1MD DAS(11)= 50.07744430D0+ 12235.1106862167D0*DTM1 2 +0.0519078250D0*DTM2 3 -0.0000298556D0*DTM3 4 -0.0000972333D0*DTM4 5 -0.0000045278D0*DTM5 6 +0.0000002861D0*DTM6 DASP(11)= (+12235.1106862167D0 2 +2.D0*0.0519078250D0*DTM1 3 -3.D0*0.0000298556D0*DTM2 4 -4.D0*0.0000972333D0*DTM3 5 -5.D0*0.0000045278D0*DTM4 6 +6.D0*0.0000002861D0*DTM5)*D1MD DO 3110 I=1,11 DAS(I)=DMOD(DAS(I),360.0D0) IF(DAS(I).LT.0.D0) DAS(I)=DAS(I)+360.0D0 3110 CONTINUE IF(IPRINT.EQ.0) RETURN C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Print astronomical elements: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! WRITE(IUN16,17004) (DAS(K),DASP(K),K=1,11) C 5000 CONTINUE WRITE(IUN16,17030) RETURN C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Format statements: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 17001 FORMAT(//6x,'Routine ETASTN, version 1996.05.25.'// 1 6x,'Astronomic elements for initial epoch '/ 2 6x,'Modified Julian date (TDT) : ',F15.4/) 17004 FORMAT(// 1 6x,'local Moontime F01',F20.11,' deg F01.',F18.11,' deg/h'/ 2 6x,'lunar longitude F02',F20.11,' deg F02.',F18.11,' deg/h'/ 3 6x,'solar longitude F03',F20.11,' deg F03.',F18.11,' deg/h'/ 4 6x,'lunar perigee F04',F20.11,' deg F04.',F18.11,' deg/h'/ 5 6x,'lunar node longit. F05',F20.11,' deg F05.',F18.11,' deg/h'/ 6 6x,'solar perigee F06',F20.11,' deg F06.',F18.11,' deg/h'/ 7 6x,'longitude Mercury F07',F20.11,' deg F07.',F18.11,' deg/h'/ 8 6x,'longitude Venus F08',F20.11,' deg F08.',F18.11,' deg/h'/ 9 6x,'longitude Mars F09',F20.11,' deg F09.',F18.11,' deg/h'/ . 6x,'longitude Jupiter F10',F20.11,' deg F10.',F18.11,' deg/h'/ 1 6x,'longitude Saturn F11',F20.11,' deg F11.',F18.11,' deg/h'/ 2) 17030 FORMAT(/6x,'***** Routine ETASTN finished the execution.'/) END C SUBROUTINE ETDDTA(IUN16,IUN27,IPRINT) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C ! C Routine ETDDTA, version 1996.05.29 Fortran 90. ! C ! C The routine ETDDTA reads a table of DDT = ET -UTC or TDT - UTC ! C from file etddt.dat. The file will be opened and after use ! C closed by the routine. ! C ! C The table on file etddt.dat has to be extended, when new data ! C are available. ! C ! C Input parameter description: ! C ---------------------------- ! C ! C IUN16: Unit number of formatted printer unit. ! C IUN27: Unit number of formmated unit, on which the table ! C of DDT has to be stored before the call of routine ! C ETDDTA. This unit will be opened by routine ETDDTA ! C as /home/hwz/eterna34/commdat/etddt.dat ! C IPRINT: Printout parameter. For IPRINT=0, nothing will be ! C written to IUN16. ! C ! C COMMON /DDT/: ! C ------------- ! C DDTTAB: Array (1..3,1..100) containing the table of year, ! C Julian date and DDT. ! C NDDTAB: Number of defined entries in table DDTTAB. ! C ! C Routine creation: 1995.12.20 by Hans-Georg Wenzel, ! C Black Forest Observatory, ! C Universitaet Karlsruhe, ! C Englerstr. 7, ! C D-76128 KARLSRUHE, ! C Germany. ! C Tel.: 0721-6082307, ! C FAX: 0721-694552. ! C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de ! C Last modification: 1996.05.29 by Hans-Georg Wenzel, ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IMPLICIT DOUBLE PRECISION (D) IMPLICIT INTEGER (I-N) CHARACTER*10 CTEXT(8),CENDT DOUBLE PRECISION DDTTAB(3,300) COMMON /DDT/ DDTTAB,NDDTAB SAVE DATA CENDT/'C*********'/ C HP-UNIX: OPEN(UNIT=27,FILE='../commdat/etddt.dat',STATUS='OLD') C MS-DOS: OPEN(UNIT=IUN27,FILE='/home/hwz/eterna34/commdat/etddt.dat', 1 STATUS='OLD') 100 READ(IUN27,17001) (CTEXT(I),I=1,8) IF(IPRINT.GT.0) WRITE(IUN16,17002) (CTEXT(I),I=1,8) IF(CTEXT(1).NE.CENDT) GOTO 100 NDDTAB=1 200 READ(IUN27,17003,END=1000) DDTTAB(1,NDDTAB),DDTTAB(2,NDDTAB), 1 DDTTAB(3,NDDTAB) IF(IPRINT.NE.0) THEN WRITE(IUN16,17004) DDTTAB(1,NDDTAB),DDTTAB(2,NDDTAB), 1 DDTTAB(3,NDDTAB) ENDIF NDDTAB=NDDTAB+1 GOTO 200 1000 NDDTAB=NDDTAB-1 CLOSE(IUN27) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Format statements: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 17001 FORMAT(8A10) 17002 FORMAT(1X,7A10,A8) 17003 FORMAT(F15.5,F15.6,F15.3) 17004 FORMAT(F15.5,F15.6,F15.3) RETURN END C SUBROUTINE ETDDTB(IUN16,IPRINT,DTUJD,DDT) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C ! C Routine ETDDTB, version 1996.05.25 Fortran 90. ! C ! C All variables with D as first character are DOUBLE PRECISION. ! C ! C Input parameter description: ! C ---------------------------- ! C ! C IUN16: Formatted printer unit. ! C IPRINT: Printout parameter. For IPRINT=0, nothing will be ! C written on unit IUN16. ! C DTUJD: Julian date of epoch (Universal time). ! C ! C Output parameter description: ! C ----------------------------- ! C ! C DDT: Difference ET - UTC resp. TDT - UTC in seconds ! C from 1955.5 until now. For epochs less 1955.5, DDT ! C is set to 31.59 s. ! C For epochs exceeding the last tabulated epoch, DDT ! C is set to the last tabulated DDT. ! C ET is Ephemeris Time. ! C TDT is Terrestrial Dynamical Time. ! C UTC is Universal Time Coordinated, as broadcasted ! C by radio or GPS satellites. ! C ! C COMMON /DDT/: ! C ------------- ! C ! C DDTTAB: Array (1..3,1..300) containing the table of year, ! C Julian date and DDT. ! C NDDTAB: Number of defined entries in table DDTTAB. ! C ! C Execution time: ! C --------------- ! C ! C 1.38 microsec per call on a 100 MHz Pentium using Lahey LF90 ! C compiler. ! C ! C Routine creation: 1995.12.20 by Hans-Georg Wenzel, ! C Black Forest Observatory, ! C Universitaet Karlsruhe, ! C Englerstr. 7, ! C D-76128 KARLSRUHE, ! C Germany. ! C Tel.: 0721-6082307, ! C FAX: 0721-694552. ! C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de ! C Last modification: 1996.05.25 by Hans-Georg Wenzel, ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IMPLICIT DOUBLE PRECISION (D) IMPLICIT INTEGER (I-N) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C COMMON /DDT/: stored table DDTTAB of DDT = TDT - UTC: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DOUBLE PRECISION DDTTAB(3,300) COMMON /DDT/ DDTTAB,NDDTAB SAVE DATA IWARN/1/,ITAB/1/ IF(DTUJD.LT.DDTTAB(2,NDDTAB)) GOTO 100 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C DTUJD exceeds last tabulated epoch DDTTAB(2,NDDTAB). ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DDT=DDTTAB(3,NDDTAB) IF(IWARN.EQ.1) WRITE(IUN16,17003) DDTTAB(1,NDDTAB) IWARN=0 RETURN C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Look at table at position ITAB. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 100 CONTINUE IF(DTUJD.GE.DDTTAB(2,ITAB).AND.DTUJD.LT.DDTTAB(2,ITAB+1)) GOTO 230 IF(DTUJD.LT.DDTTAB(2,ITAB)) THEN ITAB=ITAB-1 IF(ITAB.GT.0) GOTO 100 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Set DDT to first tabulated value and return: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ITAB=1 DDT=DDTTAB(3,1) RETURN ENDIF IF(DTUJD.GT.DDTTAB(2,ITAB+1)) THEN ITAB=ITAB+1 IF(ITAB.LT.NDDTAB) GOTO 100 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Set DDT to last tabulated value and return: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ITAB=NDDTAB DDT=DDTTAB(3,NDDTAB) RETURN ENDIF C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Interpolate table between position ITAB and ITAB+1: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 230 DDT=(DDTTAB(3,ITAB+1)*(DTUJD-DDTTAB(2,ITAB))-DDTTAB(3,ITAB)* 1 (DTUJD-DDTTAB(2,ITAB+1)))/(DDTTAB(2,ITAB+1)-DDTTAB(2,ITAB)) RETURN C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Format statements: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 17001 FORMAT(//' Routine ETDDTB.FOR, version 1996.05.25.'// 1' List of tables:'// 2' No. Juld DTX DTY'//) 17002 FORMAT(I10,2F15.5,F10.3) 17003 FORMAT(/ 1' ***** Warning from routine ETDDTB.FOR, version 1996.05.25.'/ 2' ***** Epoch exceeds the last tabulated value:',F10.5/ 3' ***** DDT of last tabulated epoch is used.'/ 4' ***** Please try to update tabels in file etddt.dat.'/) END C SUBROUTINE PREDIN(IUN15,IUN16,IPRINT) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C ! C Routine PREDIN, version 1996.05.25 Fortran 90. ! C ! C The routine reads the control parameter file *.INI and returns ! C the control parameters via COMMON /CONTROL3/ and /CONTROL4/. ! C ! C Description of input parameters: ! C -------------------------------- ! C ! C IUN15: Unit number of formatted control parameter file. ! C IUN16: Unit number of formatted printer file. ! C IPRINT: Print parameter. For IPRINT = 0, nothing will be ! C written on IUN16. ! C ! C Description of COMMON block /CONTROL3/: ! C --------------------------------------- ! C ! C DLAT: Stations latitude in degree, referring to WGS84. ! C DLON: Stations longitude in degree, positiv east of ! C Greenwich, referring to WGS84. ! C DH: Ellipsoidal height of the station in meter ! C referring to WGS84. ! C DGRAV: Gravity of the station in m/s**2 (necessary for ! C tidal tilt only). ! C DAZ: Azimuth of the earth tide sensor in degree (only ! C for tilt and horizontal strain). ! C DFRA: Lowest frequency within wave group in deg/h. ! C DFRE: Highest frequency within wave group in deg/h. ! C DG0: amplitude factor. ! C DPHI0: phase lead in degree. ! C DATLIM: Threshold for data snooping. ! C DAMIN: Amplitude threshold for the tidal potential ! C catalogue in m**2/s**2. ! C DPOLTC: Pole tide amplitude factor. ! C DLODTC: LOD tide amplitude factor. ! C IDTSEC: Sampling interval in seconds. ! C ! C Description of COMMON block /CONTROL4/: ! C --------------------------------------- ! C ! C IC: Earth tide component. ! C IR: ! C ITY: Year of the initial epoch. ! C ITM: Month of the initial epoch. ! C ITD: Day of the initial epoch. ! C ITH: Hour (UTC) of the initial epoch. ! C IDA: ! C KFILT: ! C IPROBS: ! C IPRLF: ! C IMODEL: Tidal potential catalogue. ! C IRIGID: ! C IHANN: ! C IQUICK: ! C ISPANH: ! C NGR: Number of wave groups. ! C NF: Number of meteorological parameters. ! C IREG: ! C CFY1: ! C CFY2: ! C CINST: Earth tide sensor name (CHARACTER*10). ! C CNSY: ! C CHEAD: ! C ! C Program creation: 1994.11.01 by Hans-Georg Wenzel, ! C Black Forest Observatory, ! C Universitaet Karlsruhe, ! C Englerstr. 7, ! C D-76128 KARLSRUHE 1. ! C Tel.: 0721-6082301. ! C FAX: 0721-694552. ! C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de ! C Last Modification: 1996.05.25 by Hans-Georg Wenzel. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IMPLICIT DOUBLE PRECISION (D) IMPLICIT INTEGER (I-N) CHARACTER CINPUT*75,CINTERN*50,CONTROL*10,CREST*64,CINST*10 CHARACTER CHEAD(10)*64 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C The following dimension statement is concerning the number of ! C meteorological parameters, which is 8 in the current program ! C version. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! PARAMETER (MAXNF=8) INTEGER IREG(MAXNF) CHARACTER CFY1(MAXNF)*10,CFY2(MAXNF)*10 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C The following dimension statements are concerning the number of ! C wavegroups to be used, which is 85 in the current program ! C version (parameter MAXWG). ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! PARAMETER (MAXWG=85) DIMENSION DFRA(MAXWG),DFRE(MAXWG),DG0(MAXWG),DPHI0(MAXWG) CHARACTER CNSY(MAXWG)*4 COMMON /CONTROL3/ DLAT,DLON,DH,DGRAV,DAZ,DFRA,DFRE,DG0,DPHI0, 1 DATLIM,DAMIN,DPOLTC,DLODTC,IDTSEC COMMON /CONTROL4/ IC,IR,ITY,ITM,ITD,ITH,IDA,KFILT,IPROBS, 1 IPRLF,IMODEL,IRIGID,IHANN,IQUICK,ISPANH,NGR,NF, 2 IREG,CFY1,CFY2,CINST,CNSY,CHEAD SAVE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Define default parameters: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IH=1 IGR=0 IF=0 DAMIN=1.D-8 DPOLTC=1.16D0 DLODTC=1.16D0 ITH=0 IMODEL=7 IRIGID=0 REWIND IUN15 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Read control record: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 100 CONTINUE READ(IUN15,17001,END=5000) CINPUT II=INDEX(CINPUT,'=') IF(II.EQ.0) GOTO 100 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Input record contains an equal sign at position II: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CONTROL=CINPUT(1:II-1) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Search for # in the same record: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! NLE=LEN(CINPUT) INBL=NLE DO 200 I=II+1,NLE IF(CINPUT(I:I).NE.'#') GOTO 200 INBL=I GOTO 210 200 CONTINUE 210 CREST=CINPUT(II+1:INBL-1) C WRITE(IUN16,17006) CREST C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Search for sensor name: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(CONTROL.NE.'SENSORNAME') GOTO 1300 CINST=CREST GOTO 100 1300 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Search for sampling interval: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(CONTROL.NE.'SAMPLERATE') GOTO 1400 WRITE(CINTERN,'(A15)') CREST READ(CINTERN,'(I15)') IDTSEC GOTO 100 1400 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Search for stations latitude: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(CONTROL.NE.'STATLATITU') GOTO 2400 WRITE(CINTERN,'(A15)') CREST READ(CINTERN,'(F15.4)') DLAT GOTO 100 2400 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Search for stations longitude: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(CONTROL.NE.'STATLONITU') GOTO 2500 WRITE(CINTERN,'(A15)') CREST READ(CINTERN,'(F15.4)') DLON GOTO 100 2500 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Search for stations height: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(CONTROL.NE.'STATELEVAT') GOTO 2600 WRITE(CINTERN,'(A15)') CREST READ(CINTERN,'(F15.4)') DH GOTO 100 2600 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Search for stations gravity: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(CONTROL.NE.'STATGRAVIT') GOTO 2700 WRITE(CINTERN,'(A15)') CREST READ(CINTERN,'(F15.4)') DGRAV GOTO 100 2700 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Search for stations azimuth: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(CONTROL.NE.'STATAZIMUT') GOTO 2800 WRITE(CINTERN,'(A15)') CREST READ(CINTERN,'(F15.4)') DAZ GOTO 100 2800 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Search for initial epoch: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(CONTROL.NE.'INITIALEPO') GOTO 2900 WRITE(CINTERN,'(A15)') CREST READ(CINTERN,'(3I5)') ITY,ITM,ITD GOTO 100 2900 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Search for tidal component: ! C ! C IC... Earth tide component to be computed. ! C IC=-1: tidal potential, geodetic coefficients ! C in m**2/s**2. ! C IC= 0: vertical tidal acceleration (gravity tide), ! C geodetic coefficients in nm/s**2 (positive ! C down). ! C IC= 1: horizontal tidal acceleration (tidal tilt) ! C in azimuth DAZ, geodetic coefficients in ! C mas = arc sec/1000. ! C IC= 2: vertical tidal displacement, geodetic ! C coefficients in mm. ! C IC= 3: horizontal tidal displacement in azimuth ! C DAZ, geodetic coefficients in mm. ! C IC= 4: vertical tidal strain, geodetic coefficients ! C in 10**-9 = nstr. ! C IC= 5: horizontal tidal strain in azimuth DAZ, ! C geodetic coefficients in 10**-9 = nstr. ! C IC= 6: areal tidal strain, geodetic coefficients ! C in 10**-9 = nstr. ! C IC= 7: shear tidal strain, geodetic coefficients ! C in 10**-9 = nstr. ! C IC= 8: volume tidal strain, geodetic coefficients ! C in 10**-9 = nstr. ! C IC= 9: ocean tides, geodetic coefficients in ! C millimeter. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(CONTROL.NE.'TIDALCOMPO') GOTO 3000 WRITE(CINTERN,'(A15)') CREST READ(CINTERN,'(I15)') IC GOTO 100 3000 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Search for tidal potential catalogue: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(CONTROL.NE.'TIDALPOTEN') GOTO 3100 WRITE(CINTERN,'(A15)') CREST READ(CINTERN,'(I15)') IMODEL GOTO 100 3100 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Search for truncation parameter: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(CONTROL.NE.'AMTRUNCATE') GOTO 3150 WRITE(CINTERN,'(A15)') CREST READ(CINTERN,'(F15.10)') DAMIN GOTO 100 3150 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Search for print parameter of tidal component development: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(CONTROL.NE.'PRINTDEVEL') GOTO 3200 WRITE(CINTERN,'(A15)') CREST READ(CINTERN,'(I15)') IR GOTO 100 3200 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Search for textheader: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(CONTROL.NE.'TEXTHEADER') GOTO 3400 IF(IH.GT.10) GOTO 3300 CHEAD(IH)=CREST IH=IH+1 GOTO 100 3300 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Search for data error search threshold: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(CONTROL.NE.'SEARDATLIM') GOTO 3400 WRITE(CINTERN,'(A15)') CREST READ(CINTERN,'(F15.4)') DATLIM IDA=1 IF(DATLIM.LE.0.D0) IDA=0 GOTO 100 3400 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Search for numerical lowpass filter to be selected: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(CONTROL.NE.'NUMHIGPASS') GOTO 3500 WRITE(CINTERN,'(A15)') CREST READ(CINTERN,'(I15)') KFILT GOTO 100 3500 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Search for print parameter for observations: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(CONTROL.NE.'PRINTOBSER') GOTO 3600 WRITE(CINTERN,'(A15)') CREST READ(CINTERN,'(I15)') IPROBS GOTO 100 3600 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Search for print parameter for observations: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(CONTROL.NE.'PRINTLFOBS') GOTO 3700 WRITE(CINTERN,'(A15)') CREST READ(CINTERN,'(I15)') IPRLF GOTO 100 3700 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Search for rigid earth model parameter: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(CONTROL.NE.'RIGIDEARTH') GOTO 3800 WRITE(CINTERN,'(A15)') CREST READ(CINTERN,'(I15)') IRIGID GOTO 100 3800 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Search for Hann-window parameter: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(CONTROL.NE.'HANNWINDOW') GOTO 3900 WRITE(CINTERN,'(A15)') CREST READ(CINTERN,'(I15)') IHANN IF(IHANN.GT.1) IHANN=1 GOTO 100 3900 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Search for quick look adjustment parameter: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(CONTROL.NE.'QUICKLOOKA') GOTO 4000 WRITE(CINTERN,'(A15)') CREST READ(CINTERN,'(I15)') IQUICK GOTO 100 4000 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Search for pole tide correction parameter: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(CONTROL.NE.'POLTIDECOR') GOTO 4100 WRITE(CINTERN,'(A15)') CREST READ(CINTERN,'(F15.5)') DPOLTC GOTO 100 4100 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Search for length of day tide correction parameter: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(CONTROL.NE.'LODTIDECOR') GOTO 4200 WRITE(CINTERN,'(A15)') CREST READ(CINTERN,'(F15.5)') DLODTC GOTO 100 4200 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Search for prediction span: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(CONTROL.NE.'PREDICSPAN') GOTO 4300 WRITE(CINTERN,'(A15)') CREST READ(CINTERN,'(I15)') ISPANH GOTO 100 4300 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Search for tidal parameters: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(CONTROL.NE.'TIDALPARAM') GOTO 4400 IGR=IGR+1 WRITE(CINTERN,'(A45)') CREST C READ(CINTERN,'(4F10.4,1X,A4)') DFRA(IGR),DFRE(IGR),DG0(IGR), READ(CINTERN,'(2F10.4,2F9.3,A5)') DFRA(IGR),DFRE(IGR),DG0(IGR), 1 DPHI0(IGR),CNSY(IGR) GOTO 100 4400 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Unknown control parameter name: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C WRITE(IUN16,17005) CONTROL GOTO 100 5000 CONTINUE NGR=IGR NF=IF IF(IPRINT.EQ.0) RETURN C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Print first part of control parameters: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! NH=IH DO 5010 IH=1,NH 5010 WRITE(IUN16,17114) CHEAD(IH) WRITE(IUN16,17111) IDTSEC,DLAT,DLON,DH,DGRAV,DAZ,DAMIN C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Print second part of control parameters: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! WRITE(IUN16,17112) IC,IR,ITY,ITM,ITD,ITH,IMODEL DO 5020 IGR=1,NGR 5020 WRITE(IUN16,17113) DFRA(IGR),DFRE(IGR),DG0(IGR),DPHI0(IGR), 1 CNSY(IGR) RETURN C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Format statements: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 17001 FORMAT(A75) 17002 FORMAT(1X,A75) 17003 FORMAT(' control parameter is: ',A10) 17004 FORMAT(' numerical parameter is: ',A10) 17005 FORMAT(' *** unknown control parameter:',A10) 17006 FORMAT(' CREST=',A15) 17111 FORMAT(/ 1 6x,'sample interval :',I10,' s'/ 2 6x,'stations latitude in degree :',F10.4/ 3 6x,'stations longitude in degree :',F10.4/ 4 6x,'stations height in meter :',F10.3/ 5 6x,'stations gravity in m/s**2 :',F10.4/ 6 6x,'stations azimuth from north in degree :',F10.4/ 7 6x,'truncation of tidal waves at (m**2/s**2) :',D10.3) 17112 FORMAT( 1 6x,'earth tide component : ',I10/ 2 6x,'print tidal component development (1=yes) : ',I10/ 3 6x,'initial epoch for tidal development : ',I4,I3,I3,I3/ 4 6x,'tidal potential catalogue : ',I10/) 17113 FORMAT(6x,'wave group : ',2F10.6,2F10.4,1X,A4) 17114 FORMAT(1X,A64) END C SUBROUTINE ETGCON(IUN16,IPRINT,DLAT,DLON,DH,DGRAV,DAZ,IC,DGK,DPK) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C ! C Routine ETGCON, version 1997.03.03 Fortran 90. ! C corrected 2005.06.28 (B.Ducarme) ! C The routine ETGCON computes the geodetic coefficients for ! C the tidal potential developments, Hartmann and Wenzel ! C normalization. ! C ! C All variables with D as first character are DOUBLE PRECISION. ! C ! C Input parameter description: ! C ---------------------------- ! C ! C IUN16: formatted line printer unit. ! C DLAT: ellipsoidal latitude in degree, referring to ! C geodetic reference system GRS80. ! C DLON: ellipsoidal longitude in degree, referring to ! C geodetic reference system GRS80, positiv east of ! C Greenwhich. ! C DH: ellipsoidal height in meter, referring to geodetic ! C reference system GRS80. ! C DGRAV: gravity in m/s**2. If DGRAV less than 9.50 m/s**2, ! C DGRAV will be overwritten by normal gravity ! C referring to geodetic reference system 1980. ! C DAZ: azimuth in degree from north direction counted ! C clockwise (necessary for tidal tilt only). ! C IC: Earth tide component to be computed. ! C IC=-1: tidal potential, geodetic coefficients ! C in m**2/s**2. ! C IC= 0: vertical tidal acceleration (gravity tide), ! C geodetic coefficients in nm/s**2 (positive ! C down). ! C IC= 1: horizontal tidal acceleration (tidal tilt) ! C in azimuth DAZ, geodetic coefficients in ! C mas = arc sec/1000. ! C IC= 2: vertical tidal displacement, geodetic ! C coefficients in mm. ! C IC= 3: horizontal tidal displacement in azimuth ! C DAZ, geodetic coefficients in mm. ! C IC= 4: vertical tidal strain, geodetic coefficients ! C in 10**-9 = nstr. ! C IC= 5: horizontal tidal strain in azimuth DAZ, ! C geodetic coefficients in 10**-9 = nstr. ! C IC= 6: areal tidal strain, geodetic coefficients ! C in 10**-9 = nstr. ! C IC= 7: shear tidal strain, geodetic coefficients ! C in 10**-9 = nstr. ! C IC= 8: volume tidal strain, geodetic coefficients ! C in 10**-9 = nstr. ! C IC= 9: ocean tides, geodetic coefficients in ! C millimeter. ! C ! C Output parameter description: ! C ----------------------------- ! C ! C DGK: array (1...25) of geodetic coefficients. ! C The geodetic coefficient of degree L and order M ! C is stored in DGK(J) with J=L*(L+1)/2+M-2. ! C DPK: array (1...25) of phases in degree. ! C The phase for degree L and order M is stored in ! C DPK(J) with J=L*(L+1)/2+M-2. ! C ! C Used routines: ! C -------------- ! C ! C ETLOVE: computes latitude dependent elastic parameters. ! C ETLEGN: computes fully normalized Legendre functions and their ! C derivatives. ! C ! C Numerical accuracy: ! C ------------------- ! C ! C The routine has been tested under operation system MS-DOS and ! C UNIX in double precision (8 byte words = 15 digits) using ! C different compilers. ! C ! C References: ! C ! C Wilhelm, H. and W. Zuern (1984): Tidal forcing field. ! C In: Landolt-Boernstein, Zahlenwerte und Funktionen aus ! C Naturwissenschaften und Technik, New series, group V, ! C Vol. 2, Geophysics of the Solid Earth, the Moon and the ! C Planets, Berlin 1984. ! C ! C Zuern, W. and H. Wilhelm (1984): Tides of the solid Earth. ! C In: Landolt-Boernstein, Zahlenwerte und Funktionen aus ! C Naturwissenschaften und Technik, New series, group V, Vol. ! C 2, Geophysics of the Solid Earth, the Moon and the Planets,! C Berlin 1984. ! C ! C Routine creation: 1988.01.29 by Hans-Georg Wenzel, ! C Black Forest Observatory, ! C Universitaet Karlsruhe, ! C Englerstr. 7, ! C D-76128 KARLSRUHE, ! C Germany. ! C Tel.: 0721-6082307, ! C FAX: 0721-694552. ! C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de ! C ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IMPLICIT DOUBLE PRECISION (D) IMPLICIT INTEGER (I-N) CHARACTER CUNIT(11)*8 DOUBLE PRECISION DGK(25),DPK(25),DGX(25),DGY(25),DGZ(25) DOUBLE PRECISION DP0(25),DP1(25) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C COMMON /LOVE/ contains gravimeter factors, LOVE-numbers, SHIDA- ! C numbers and tilt factors for degree 2...4 at latitude DLAT: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C COMMON /CONST/: ! C DPI... 3.1415.... ! C DPI2... 2.D0*DPI ! C DRAD... DPI/180.D0 ! C DRO... 180.D0/DPI ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! COMMON /CONST/ DPI,DPI2,DRAD,DRO COMMON /UNITS/ CUNIT,IC2 SAVE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Definition of parameters of Geodetic Reference System 1980. ! C DEA is major semi axis in meter. ! C DEE is square of first excentricity (without dimension). ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DATA DEA/6378136.3D0/,DEE/6.69439795140D-3/ IF(IPRINT.GT.0) WRITE(IUN16,17000) DEA,DEE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C DCLAT is cos and DSLAT is sin of ellipsoidal latitude. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DCLAT=DCOS(DLAT*DRAD) DSLAT=DSIN(DLAT*DRAD) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Compute normal gravity in m/s**2: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Compute ellipsoidal curvature radius DN in meter. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DN=DEA/DSQRT(1.D0-DEE*DSLAT**2) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Compute geocentric latitude DPSI in degree: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DPSI=DRO*DATAN(((DN*(1.D0-DEE)+DH)*DSLAT)/((DN+DH)*DCLAT)) DTHET=90.D0-DPSI DCT=DCOS(DTHET*DRAD) DST=DSIN(DTHET*DRAD) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Compute fully normalized spherical harmonics: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CALL ETLEGN(DCT,DST,LMAX,DP0,DP1) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Compute geocentric radius DR in meter: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DR=DSQRT((DN+DH)**2*DCLAT**2+(DN*(1.D0-DEE)+DH)**2*DSLAT**2) IF(IPRINT.GT.0) WRITE(IUN16,17001) DLAT,DPSI,DLON,DH,DGRAV,DR,IC, 1 DAZ C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C DGRAV*10.**9:nm/s**2,DRO*3600.*10.**3:radian to mas ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DF=DRO*3.600D-3/DGRAV C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Compute latitude dependent elastic parameters from Wahr-Dehant- ! C Zschau model: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CALL ETLOVE(IUN16,IPRINT,DLAT,DH) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C DCPSI is cos and DSPSI is sin of geocentric latitude. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DCPSI=DCOS(DPSI*DRAD) DSPSI=DSIN(DPSI*DRAD) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Compute spherical geodetic coefficients. ! C DGK contains coefficients for potential in m**2/s**2! C DGX contains coefficients for north accelerations in nm/s**2. ! C DGY contains coefficients for east accelerations in nm/s**2. ! C DGZ contains coefficients for vertical accelerations in nm/s**2. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DRDA=DR/DEA DO 10 LI=2,LMAX DRDADL=DRDA**LI DO 10 MI=0,LI J=LI*(LI+1)/2+MI-2 DGK(J)= DRDADL*DP0(J) DGX(J)=-1.D0*DRDADL/DR*DP1(J)*1.D9 DGY(J)= DRDADL*DBLE(MI)/(DR*DST)*DP0(J)*1.D9 DGZ(J)= DRDADL*DBLE(LI)/DR*DP0(J)*1.D9 10 CONTINUE DO 20 I=1,25 20 DPK(I)=0.D0 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Compute geodetic coefficients for tidal acceleration vector ! C orientated to ellipsoidal coordinate system stored in ! C DGX (north), DGY (east) and DGZ (upwards), all in nm/s**2. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DCDLAT=DCLAT*DCPSI+DSLAT*DSPSI DSDLAT=DSLAT*DCPSI-DCLAT*DSPSI DO 50 I=1,25 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 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C IC=-1, compute geodetic coefficients for tidal potential ! C (m**2/s**2). ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 100 CONTINUE GOTO 2000 200 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C IC=0, compute geodetic coefficients for vertical component ! C (gravity tide in nm/s**2). ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DO 210 I=1,25 DGK(I)=DGZ(I) 210 DPK(I)=180.0D0 GOTO 2000 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C IC=1, compute geodetic coefficients for horizontal component ! C (tidal tilt) in azimuth DAZ, in mas. ! C DF:mas/(nm/s**2), DGX(I),DGY(I): nm/s**2 ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C IC=2, compute geodetic coefficients for vertical displacement ! C in mm. ! C DGK(I):m**2/s**2, DGRAV:m/s**2, 10**3 conversion to mm ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 400 CONTINUE DFAK=1.D3/DGRAV DO 410 I=1,12 DGK(I)=DGK(I)*DHLAT(I)*DFAK 410 DPK(I)=0.0D0 WRITE(IUN16,*) '*****The component',IC,' has never been tested !' GOTO 2000 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C IC=3, compute geodetic coefficients for horizontal displacement ! C in azimuth DAZ in mm. ! C DGRAV*10.**9:nm/s**2,10.**3:conversion to mm (corr. 2004.02.18) ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 500 CONTINUE DFAK=1.D-6*DR/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(IUN16,*) '*****The component',IC,' has never been tested !' GOTO 2000 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C IC=4, compute geodetic coefficients for vertical strain at the ! C Earth's deformed surface in 10**-9 units = nstr. ! C We use a spherical approximation for the vertical strain, ! C i.e. eps(rr) , and a POISSON ratio of 0.25 (see ZUERN and ! C WILHELM 1984, p. 282). ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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*DR) DO 620 I=4,7 620 DGK(I)=DGK(I)*DFAK*(2.D0*DHLAT(I)-3.D0*4.D0*DLLAT(I))/(DGRAV*DR) DO 630 I=8,12 630 DGK(I)=DGK(I)*DFAK*(2.D0*DHLAT(I)-4.D0*5.D0*DLLAT(I))/(DGRAV*DR) DO 640 I=1,12 640 DPK(I)=0.0D0 GOTO 2000 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C IC=5, compute geodetic coefficients for horizontal strain ! C in azimuth DAZ, in 10**-9 units. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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/(DR*DGRAV) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Real part is stored in DGX, imaginary part is stored in DGY. ! C Formulas were given by Dr. W. Zuern, BFO Schiltach (personal ! C communication) and tested against horizontal strain computed ! C (with lower precision) by program ETIDEL (made by Bilham). ! C Results agreed to 0.3 % and 0.1 degree for most of the waves, ! C except for 2N2 and L2 (deviation of 3 %). ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C IC=6, compute geodetic coefficients for areal strain ! C in 10**-9 units = nstr. ! C We use a spherical approximation for the aereal strain, ! C i.e. eps(t,t) + eps(l,l), (see ZUERN and WILHELM 1984, ! C p. 282). ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 800 CONTINUE DO 810 I=1,3 810 DGK(I)=DGK(I)*(2.D0*DHLAT(I)-2.D0*3.D0*DLLAT(I))/(DGRAV*DR)*1.D9 DO 820 I=4,7 820 DGK(I)=DGK(I)*(2.D0*DHLAT(I)-3.D0*4.D0*DLLAT(I))/(DGRAV*DR)*1.D9 DO 830 I=8,12 830 DGK(I)=DGK(I)*(2.D0*DHLAT(I)-4.D0*5.D0*DLLAT(I))/(DGRAV*DR)*1.D9 DO 840 I=1,12 840 DPK(I)=0.0D0 GOTO 2000 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C IC=7, compute geodetic coefficients for shear tidal strain ! C at the Earth's deformed surface in 10**-9 units = nstr. ! C We use a spherical approximation, i.e. eps(t,l) ! C Attention: this component has never been tested !!!! ! C corr. 2005/06/29 ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 900 CONTINUE DTHETA=(90.D0-DPSI)*DRAD c DAZR=(DAZ+180.D0)*DRAD c DCAZ =DCOS(DAZR) c DSAZ =DSIN(DAZR) c DSAZ2=DSIN(2.D0*DAZR) c 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/(DR*DGRAV) DGY(1)=0.D0 c DGY(2)=2.D0*DLLAT(2)*(2.D0*DCOTT2-DCOTT)*DCSTS/DST DGY(2)=2.D0*DLLAT(2)*(2.D0*DCOTT2-DCOTT)/DST c DGY(3)=4.D0*DLLAT(3)*DCOTT*DCSTS/DST DGY(3)=4.D0*DLLAT(3)*DCOTT/DST DGY(4)=0.D0 c DGY(5)=-20.D0*DLLAT(5)*DCT*DCSTS/(5.D0*DCT2-1.D0) DGY(5)=-20.D0*DLLAT(5)*DCT/(5.D0*DCT2-1.D0) c DGY(6)=-4.D0*DLLAT(6)*(DCOTT-1.D0/DCOTT)*DCSTS/DST DGY(6)=-4.D0*DLLAT(6)*(DCOTT-1.D0/DCOTT)/DST c DGY(7)=12.D0*DLLAT(7)*DCOTT*DCSTS/DST DGY(7)=12.D0*DLLAT(7)*DCOTT/DST DGY(8)=0.D0 c DGY(9)=DLLAT(9)*3.D0/DCT*(1.D0+2.D0/(7.D0*DCT2-3.D0))*DSAZ2 DGY(9)=DLLAT(9)*3.D0/DCT*(1.D0+2.D0/(7.D0*DCT2-3.D0)) c DGY(10)=-DLLAT(10)*6.D0*DCT/DST**2*(1.D0-4.D0/(7.D0*DCT2-1.D0))* c 1 DSAZ2 DGY(10)=-DLLAT(10)*6.D0*DCT/DST**2*(1.D0-4.D0/(7.D0*DCT2-1.D0)) c DGY(11)=DLLAT(11)*3.D0/DCT*(3.D0-2.D0/DST2)*DSAZ2 DGY(11)=DLLAT(11)*3.D0/DCT*(3.D0-2.D0/DST2) c DGY(12)=DLLAT(12)*12.D0*DCT/DST2*DSAZ2 DGY(12)=DLLAT(12)*12.D0*DCT/DST2 DO 910 I=1,12 DGK(I)=DGK(I)*DGY(I)*DFAK 910 DPK(I)=0.D0 WRITE(IUN16,*) ' ***** The shear strain has never been tested !' GOTO 2000 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C IC=8, compute geodetic coefficients for volume strain ! C at the Earth's deformed surface in 10**-9 units = nstr. ! C We use a spherical approximation, i.e. eps(t,t)+eps(l,l). ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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*DR) DO 1020 I=4,7 1020 DGK(I)=DGK(I)*DFAK*(2.D0*DHLAT(I)-3.D0*4.D0*DLLAT(I))/(DGRAV*DR) DO 1030 I=8,12 1030 DGK(I)=DGK(I)*DFAK*(2.D0*DHLAT(I)-4.D0*5.D0*DLLAT(I))/(DGRAV*DR) DO 1040 I=1,12 1040 DPK(I)=0.0D0 GOTO 2000 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C IC=9, compute geodetic coefficients for static ocean tides in mm.! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1100 CONTINUE DFAK=1.D3/DGRAV DO 1110 I=1,25 DGK(I)=DGK(I)*DFAK 1110 DPK(I)=0.0D0 2000 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Print geodetic coefficients: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(IPRINT.EQ.0) RETURN WRITE(IUN16,17003) IC,DAZ,(DGK(I),CUNIT(IC2),DPK(I),I=1,12) if(ic.gt.0.and.ic.lt.9)go to 5000 WRITE(IUN16,17004) (DGK(I),CUNIT(IC2),DPK(I),I=13,25) 5000 CONTINUE IF(IPRINT.GT.0) WRITE(IUN16,17005) RETURN C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Format statements: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 17000 FORMAT(' Routine ETGCON, version 1997.03.03.'// 1' Computation of geodetic coefficients'// 3' Parameters of Geodetic Reference System 1980:'/ 4' Major semi axis ',F12.0,' m'/ 5' 1. excentricity ',F12.8/) 17001 FORMAT(' Station parameters:'// 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'//) 17003 FORMAT(/' Geodetic coefficients and phases for component',I4/ 1' azimuth:',F12.6,' degree'// 2' GC 2,0',F14.8,2X,A8,2X,F14.6,' deg'/ 3' GC 2,1',F14.8,2X,A8,2X,F14.6,' deg'/ 4' GC 2,2',F14.8,2X,A8,2X,F14.6,' deg'/ 5' GC 3,0',F14.8,2X,A8,2X,F14.6,' deg'/ 6' GC 3,1',F14.8,2X,A8,2X,F14.6,' deg'/ 7' GC 3,2',F14.8,2X,A8,2X,F14.6,' deg'/ 8' GC 3,3',F14.8,2X,A8,2X,F14.6,' deg'/ 9' GC 4,0',F14.8,2X,A8,2X,F14.6,' deg'/ *' GC 4,1',F14.8,2X,A8,2X,F14.6,' deg'/ 1' GC 4,2',F14.8,2X,A8,2X,F14.6,' deg'/ 2' GC 4,3',F14.8,2X,A8,2X,F14.6,' deg'/ 3' GC 4,4',F14.8,2X,A8,2X,F14.6,' deg') 17004 FORMAT( 1' GC 5,0',F14.8,2X,A8,2X,F14.6,' deg'/ 2' GC 5,1',F14.8,2X,A8,2X,F14.6,' deg'/ 3' GC 5,2',F14.8,2X,A8,2X,F14.6,' deg'/ 4' GC 5,3',F14.8,2X,A8,2X,F14.6,' deg'/ 5' GC 5,4',F14.8,2X,A8,2X,F14.6,' deg'/ 6' GC 5,5',F14.8,2X,A8,2X,F14.6,' deg'/ 7' GC 6,0',F14.8,2X,A8,2X,F14.6,' deg'/ 8' GC 6,1',F14.8,2X,A8,2X,F14.6,' deg'/ 9' GC 6,2',F14.8,2X,A8,2X,F14.6,' deg'/ *' GC 6,3',F14.8,2X,A8,2X,F14.6,' deg'/ 1' GC 6,4',F14.8,2X,A8,2X,F14.6,' deg'/ 2' GC 6,5',F14.8,2X,A8,2X,F14.6,' deg'/ 3' GC 6,6',F14.8,2X,A8,2X,F14.6,' deg'/) 17005 FORMAT(/6x,'***** Routine ETGCON finished the execution.'/) END C SUBROUTINE ETGREN(IUN16,DJULD,ITY,ITM,ITD,DTH,NERR) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C ! C Routine ETGREN, version 1996.05.25 Fortran 90. ! C ! C The routine ETGREN computes Gregorian date from given Julian ! C date. ETGREN is a modified version of routine CALDAT given in ! C Fortran by ! C Press, W.H., S.A. Teukolsky, W.T. Vetterling and B.P. Flannery ! C (1992): Numerical recipes in FORTRAN: the art of scientific ! C computing. Second edition, Cambridge University Press ! C 1992. ! C ! C The routine has been tested and found to be correct between ! C years -3000 and +3000. ! C ! C Input parameter description: ! C ---------------------------- ! C ! C DJULD: Julian date (DOUBLE PRECISION). ! C ! C Output parameter description: ! C ----------------------------- ! C ! C ITY: year (INTEGER). ! C ITM: month (INTEGER). ! C ITD: day (INTEGER). ! C DTH: hour (DOUBLE PRECISION). ! C NERR: error code, counts the number of errors which ! C happened during the execution of routine ETGREN. ! C ! C Execution time: ! C --------------- ! C ! C 2.47 microsec per call on a Pentium 100 MHz using Lahey LF90 ! C compiler. ! C ! C Routine creation: 1995.11.04 by Hans-Georg Wenzel, ! C Black Forest Observatory, ! C Universitaet Karlsruhe, ! C Englerstr. 7, ! C D-76128 KARLSRUHE, ! C Germany. ! C Tel.: 0721-6082301. ! C FAX: 0721-694552. ! C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de ! C Last modification: 1996.05.25 by Hans-Georg Wenzel, ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IMPLICIT DOUBLE PRECISION (D) IMPLICIT INTEGER (I-N) SAVE DATA DGREG/2299160.499999D0/ NERR=0 IF(DJULD.LT.625307.D0.OR.DJULD.GT.2817153.D0) THEN NERR=1 WRITE(IUN16,17050) DJULD ENDIF JULIAN=INT(DJULD) DTH=12.D0 DFRAC=DJULD-DBLE(JULIAN) DTH=DTH+DFRAC*24.D0 IF(DTH.GE.23.9999D0) THEN DTH=DTH-24.D0 JULIAN=JULIAN+1 ENDIF IF(DJULD.GT.DGREG) THEN C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Cross over from Gregorian calendar procudes this correction: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! JALPHA=INT(((JULIAN-1867216)-0.25D0)/36524.25D0) JA=JULIAN+1+JALPHA-INT(0.25D0*DBLE(JALPHA)) ELSE JA=JULIAN ENDIF JB=JA+1524 JC=INT(6680.D0+((JB-2439870)-122.1D0)/365.25D0) JD=365*JC+INT(0.25D0*DBLE(JC)) JE=INT((JB-JD)/30.6001D0) ITD=JB-JD-INT(30.6001D0*DBLE(JE)) C ITM=JE-1 IF(ITM.GT.12) ITM=ITM-12 ITY=JC-4715 IF(ITM.GT.2) ITY=ITY-1 RETURN C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Format statements: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 17050 FORMAT(/' *****Error in routine ETGREN.FOR, version 1996.05.25.'/ 1' *****Julian date is:',F20.6/ 1' *****Year is less -3000 or greater +3000.'/ 2' *****Routine ETGREN has not been tested for this case.'/ 3' *****Routine ETGREN continues the execution.'/) END C SUBROUTINE ETJULN(IUN16,ITY,ITM,ITD,DTH,DJULD) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C ! C Routine ETJULN, version 1996.05.25 Fortran 90. ! C ! C The routine ETJULN computes the Julian date and the modified ! C Julian date. ETJULN is a modified version of routine MJD given ! C in PASCAL by Montenbruck and Pfleger (see below). ! C ! C The routine is valid for every date since year -4713. ! C Comparison with reference values between years -1410 and +3200 ! C from JPL was successfully. ! C ! C Reference: Montenbruck, O. and T. Pfleger (1989): Astronomie mit ! C dem Personal Computer. Springer Verlag, Berlin 1989. ! C ! C Input parameter description: ! C ---------------------------- ! C ! C ITY: Year (INTEGER). ! C ITM: Month (INTEGER). ! C ITD: Day (INTEGER). ! C DTH: Hour (DOUBLE PRECISION). ! C ! C Output parameter description: ! C ----------------------------- ! C ! C DJULD: Julian date (DOUBLE PRECISION). ! C 16. April -1410, 0.00 H is DJULD = 1206160.5D0 ! C 31. January -1100, 0.00 H is DJULD = 1319312.5D0 ! C 24. January -0800, 0.00 H is DJULD = 1428880.5D0 ! C 17. January -0500, 0.00 H is DJULD = 1538448.5D0 ! C 10. January -0200, 0.00 H is DJULD = 1648016.5D0 ! C 03. January 100, 0.00 H is DJULD = 1757584.5D0 ! C 29. February 400, 0.00 H is DJULD = 1867216.5D0 ! C 20. December 699, 0.00 H is DJULD = 1976720.5D0 ! C 15. February 1000, 0.00 H is DJULD = 2086352.5D0 ! C 08. February 1300, 0.00 H is DJULD = 2195920.5D0 ! C 11. February 1600, 0.00 H is DJULD = 2305488.5D0 ! C 06. February 1900, 0.00 H is DJULD = 2415056.5D0 ! C 01. January 1988, 0.00 H is DJULD = 2447161.5D0 ! C 01. February 1988, 0.00 H is DJULD = 2447192.5D0 ! C 29. February 1988, 0.00 H is DJULD = 2447220.5D0 ! C 01. March 1988, 0.00 H is DJULD = 2447221.5D0 ! C 01. February 2200, 0.00 H is DJULD = 2524624.5D0 ! C 27. January 2500, 0.00 H is DJULD = 2634192.5D0 ! C 23. January 2800, 0.00 H is DJULD = 2743760.5D0 ! C 22. December 3002, 0.00 H is DJULD = 2817872.5D0 ! C ! C To obtain the modified Julian date, subtract 2400000.5 from ! C DJULD. ! C ! C Execution time: ! C --------------- ! C ! C 2.42 microsec per call on a Pentium 100 MHz using Lahex LF90 ! C compiler. ! C ! C Routine creation: 1992.09.19 by Hans-Georg Wenzel, ! C Black Forest Observatory, ! C Universitaet Karlsruhe, ! C Englerstr. 7, ! C D-76128 KARLSRUHE, ! C Germany. ! C Tel.: 0721-6082301. ! C FAX: 0721-694552. ! C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de ! C Last modification: 1996.05.25 by Hans-Georg Wenzel. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IMPLICIT DOUBLE PRECISION (D) SAVE ITYY=ITY IF(ITM.LT.1) GOTO 5000 IF(ITM.GT.12) GOTO 5010 ITMM=ITM DA=10000.0D0*ITYY+100.0D0*ITMM+ITD IF(ITMM.LE.2) THEN ITMM=ITMM+12 ITYY=ITYY-1 END IF IF(DA.LE.15821004.1D0) THEN DB=-2+(ITYY+4716)/4-1179 ELSE DB=ITYY/400-ITYY/100+ITYY/4 ENDIF DA=365.0D0*DBLE(ITYY)-679004.0D0 DJULD=DA+DB+INT(30.6001D0*(ITMM+1))+DBLE(ITD)+DTH/24.D0 1 +2400000.5D0 RETURN 5000 WRITE(IUN16,17050) ITY,ITM,ITD,DTH STOP 5010 WRITE(IUN16,17051) ITY,ITM,ITD,DTH STOP C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Format statements: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 17050 FORMAT(/' *****Error in routine ETJULN, version 1996.05.25.'/ 1' *****Month is less 1:',2X,3I4,F12.3/ 2' *****Routine ETJULN stops the execution.'/) 17051 FORMAT(/' *****Error in routine ETJULN, version 1996.05.25.'/ 1' *****Month is greater 12:',2X,3I4,F12.3/ 2' *****Routine ETJULN stops the execution.'/) END C SUBROUTINE ETLEGN(DCT,DST,LMAX,DP0,DP1) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C ! C Routine ETLEGN, version 1996.05.25 Fortran 90. ! C ! C The routine computes the fully normalized Legendre functions ! C and their derivatives complete to degree and order 6 by explicit ! C formulas. ! C ! C Input parameter description: ! C ---------------------------- ! C ! C DCT: DOUBLE PRECISION COS of polar distance theta, for ! C which the fully normalized associated Legendre ! C functions will be computed. ! C DST: DOUBLE PRECISION SIN of polar distance theta, for ! C which the fully normalized associated Legendre ! C functions will be computed. ! C ! C Output parameter desription: ! C ----------------------------- ! C ! C LMAX: Maximum degree and order, for which the fully ! C normalized associated Legendre functions will be ! C computed. LMAX is equal to 6. ! C DP0: DOUBLE PRECISION array of fully normalized Legendre ! C functions. The fully normalized Legendre function ! C of degree L and order M is stored in ! C DP0(J) WITH J=L*(L+1)/2+M+1. ! C DP1: DOUBLE PRECISION array of first derivatives of the ! C fully normalized Legendre functions to polar ! C distance theta. The first derivative of fully ! C normalized Legendre function of degree L and order ! C M is stored in DP1(J) WITH J=L*(L+1)/2+M-2. ! C ! C Example for theta = 30 degree: ! C ! C J L M DP0(L+1,M+1) DP1(L+1,M*1) ! C ! C 1 2 0 1.39754248593737 2.90473750965556 ! C 2 2 1 -1.67705098312484 1.93649167310371 ! C 3 2 2 0.48412291827593 -1.67705098312484 ! C 4 3 0 0.85923294280422 5.45686207907072 ! C 5 3 1 -2.22775461507770 0.35078038001005 ! C 6 3 2 1.10926495933118 -3.20217211436237 ! C 7 3 3 -0.26145625829190 1.35856656995526 ! C 8 4 0 0.07031250000000 7.30708934443120 ! C 9 4 1 -2.31070453947492 -3.55756236768943 ! C 10 4 2 1.78186666957014 -3.63092188706945 ! C 11 4 3 -0.67928328497763 3.13747509950278 ! C 12 4 4 0.13865811991640 -0.96065163430871 ! C 13 5 0 -0.74051002865529 7.19033890096581 ! C 14 5 1 -1.85653752113519 -8.95158333012718 ! C 15 5 2 2.29938478949397 -1.85857059805883 ! C 16 5 3 -1.24653144252643 4.78747153809058 ! C 17 5 4 0.39826512815546 -2.52932326844337 ! C 18 5 5 -0.07271293151948 0.62971245879506 ! C 19 6 0 -1.34856068213155 4.35442243247701 ! C 20 6 1 -0.95021287641141 -14.00557979016896 ! C 21 6 2 2.47470311782905 2.56294916449777 ! C 22 6 3 -1.85592870532597 5.20453026842398 ! C 23 6 4 0.81047568870385 -4.55019988574613 ! C 24 6 5 -0.22704605589841 1.83519142087945 ! C 25 6 6 0.03784100931640 -0.39325530447417 ! C ! C Execution time: ! C --------------- ! C ! C 0.00006 sec per call of ETLEGN on 80486 DX4 100MHZ with NDEG=6. ! C ! C Program creation: 1995.03.23 by Hans-Georg Wenzel, ! C Geodaetisches Institut, ! C Universitaet Karlsruhe, ! C Englerstr. 7, ! C D-76128 KARLSRUHE, ! C Germany. ! C Tel.: 0721-6082301. ! C FAX: 0721-694552. ! C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de ! C Last modification: 1996.05.25 by Hans-Georg Wenzel. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IMPLICIT DOUBLE PRECISION (D) DOUBLE PRECISION DP0(25),DP1(25) SAVE LMAX=6 DST2=DST*DST DCT2=DCT*DCT DST3=DST2*DST DCT3=DCT2*DCT DST4=DST3*DST DCT4=DCT3*DCT DST5=DST4*DST DCT5=DCT4*DCT DST6=DST5*DST DCT6=DCT5*DCT C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Compute fully normalized Legendre functions: ! C Degree 2: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DP0(01)= DSQRT(5.D0/4.D0)*(3.D0*DCT2-1.D0) DP0(02)= DSQRT(15.D0)*DCT*DST DP0(03)= DSQRT(15.D0/4.D0)*DST2 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Degree 3: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DP0(04)= DSQRT(7.D0/4.D0)*DCT*(5.D0*DCT2-3.D0) DP0(05)= DSQRT(21.D0/8.D0)*DST*(5.D0*DCT2-1.D0) DP0(06)= DSQRT(105.D0/4.D0)*DST2*DCT DP0(07)= DSQRT(35.D0/8.D0)*DST3 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Degree 4: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DP0(08)= 3.D0/8.D0*(3.D0-30.D0*DCT2+35.D0*DCT4) DP0(09)= DSQRT(45.D0/8.D0)*DST*DCT*(7.D0*DCT2-3.D0) DP0(10)= DSQRT(45.D0/16.D0)*(-1.D0+8.D0*DCT2-7.D0*DCT4) DP0(11)= DSQRT(315.D0/8.D0)*DST3*DCT DP0(12)= DSQRT(315.D0/64.D0)*DST4 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Degree 5: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DP0(13)= DSQRT(11.D0/64.D0)*DCT*(15.D0-70.D0*DCT2+63.D0*DCT4) DP0(14)= DSQRT(165.D0/64.D0)*DST*(1.D0-14.D0*DCT2+21.D0*DCT4) DP0(15)= DSQRT(1155.D0/16.D0)*DCT*(-1.D0+4.D0*DCT2-3.D0*DCT4) DP0(16)= DSQRT(385.D0/128.D0)*DST3*(9.D0*DCT2-1.D0) DP0(17)= DSQRT(3465.D0/64.D0)*DCT*DST4 DP0(18)= DSQRT(693.D0/128.D0)*DST5 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Degree 6: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DP0(19)= DSQRT(13.D0/256.D0)*(-5.D0+105.D0*DCT2-315.D0*DCT4 1 +231.D0*DCT6) DP0(20)= DSQRT(273.D0/64.D0)*DST*DCT*(5.D0-30.D0*DCT2 1 +33.D0*DCT4) DP0(21)= DSQRT(2730.D0/1024.D0)*(1.D0-19.D0*DCT2+51.D0*DCT4 1 -33.D0*DCT6) DP0(22)= DSQRT(2730.D0/256.D0)*DST3*DCT*(-3.D0+11.D0*DCT2) DP0(23)= DSQRT(819.D0/256.D0)*(-1.D0+13.D0*DCT2-23.D0*DCT4 1 +11.D0*DCT6) DP0(24)= DSQRT(18018.D0/256.D0)*DST5*DCT DP0(25)= DSQRT(6006.D0/1024.D0)*DST6 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Compute derivations with respect to theta: ! C Degree 2: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DP1(01)=-DSQRT(45.D0)*DST*DCT DP1(02)= DSQRT(15.D0)*(1.D0-2.D0*DST2) DP1(03)= DSQRT(15.D0)*DST*DCT C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Degree 3: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DP1(04)=-DSQRT(63.D0/4.D0)*DST*(5.D0*DCT2-1.D0) DP1(05)= DSQRT(21.D0/8.D0)*DCT*(4.D0-15.D0*DST2) DP1(06)=-DSQRT(105.D0/4.D0)*DST*(1.D0-3.D0*DCT2) DP1(07)= DSQRT(315.D0/8.D0)*DST2*DCT C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Degree 4: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DP1(08)=-15.D0/2.D0*(7.D0*DCT2-3.D0)*DST*DCT DP1(09)= DSQRT(45.D0/8.D0)*(3.D0-27.D0*DCT2+28.D0*DCT4) DP1(10)=-DSQRT(45.D0)*(4.D0-7.D0*DCT2)*DST*DCT DP1(11)= DSQRT(315.D0/8.D0)*DST2*(4.D0*DCT2-1.D0) DP1(12)= DSQRT(315.D0/4.D0)*DST3*DCT C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Degree 5: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DP1(13)=-DSQRT(2475.D0/64.D0)*DST*(1.D0-14.D0*DCT2+21.D0*DCT4) DP1(14)= DSQRT(165.D0/64.D0)*DCT*(29.D0-126.D0*DCT2 1 +105.D0*DCT4) DP1(15)=-DSQRT(1155.D0/16.D0)*DST*(-1.D0+12.D0*DCT2-15.D0*DCT4) DP1(16)= DSQRT(3465.D0/128.D0)*DST2*DCT*(15.D0*DCT2-7.D0) DP1(17)=-DSQRT(3465.D0/64.D0)*DST*(1.D0-6.D0*DCT2+5.D0*DCT4) DP1(18)= DSQRT(17325.D0/128.D0)*DCT*DST4 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Degree 6: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DP1(19)=-DSQRT(5733.D0/64.D0)*DST*DCT*(5.D0-30.D0*DCT2 1 +33.D0*DCT4) DP1(20)=-DSQRT(273.D0/64.D0)*(5.D0-100.D0*DCT2+285.D0*DCT4 1 -198.D0*DCT6) DP1(21)=-DSQRT(1365.D0/128.D0)*DST*DCT*(-19.D0+102.D0*DCT2 1 -99.D0*DCT4) DP1(22)= DSQRT(12285.D0/128.D0)*DST2*(1.D0-15.D0*DCT2 1 +22.D0*DCT4) DP1(23)=-DSQRT(819.D0/64.D0)*DCT*DST*(13.D0-46.D0*DCT2 1 +33.D0*DCT4) DP1(24)= DSQRT(9009.D0/128.D0)*DST4*(6.D0*DCT2-1.D0) DP1(25)= DSQRT(27027.D0/128.D0)*DST5*DCT RETURN END C SUBROUTINE ETLOVE(IUN16,IPRINT,DLAT,DELV) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C ! C Routine ETLOVE, version 1996.05.25 Fortran 90. ! C ! C The routine computes latitude dependent LOVE-numbers DH, DK, ! C SHIDA-numbers DL, gravimeter factors DG and tilt factors DT ! C using the so-called Wahr-Dehant-Zschau model. ! C ! C Body tide amplitude factors for Wahr-Dehant-Zschau model. ! C The NDFW resonance is approximated by ! C ! C G(RES) = GLAT - GR*(DOM - DOM0)/(DOMR - DOM). ! C ! C similar equations hold for the other parameters. ! C ! C Gravimetric amplitude factors, LOVE numbers h and k for degree ! C 0...3 have been taken from Dehant (1987), Table 7, 8 and 9 ! C for an elliptical, uniformly rotating, oceanless Earth with ! C liquid outer core and inelastic mantle (PREM Earth model with ! C inelastic mantle from Zschau) and for the fourth degree from ! C Dehant et. al (1989), Table 6. The resonance factors GR have ! C been computed to fit the difference between body tide amplitude ! C factors at O1 and PSI1 from Dehant (1987), PREM model with ! C elastic mantle (Table 1...3). The NDFW resonance frequency is ! C 15.073729 degree per hour = 1.004915267 CPD UT, taken from ! C Wahr (1981) (because it is not given in Dehant's papers). ! C ! C Input parameter description: ! C ---------------------------- ! C ! C IUN16: formatted line printer unit. ! C IPRINT: printout parameter. For IPRINT=1, the computed ! C Love- and Shida- number s will be printed. ! C DLAT: ellipsoidal latitude in degree. ! C DELV: ellipsoidal height in meter. ! C ! C Description of COMMON /LOVE/: ! C ----------------------------- ! C ! C DOM0: frequency of O1 in degree per hour. ! C DOMR: frequency of the FCN eigenfrequency in degree per ! C hour. ! C DGLAT: array(1..12) containing the gravimetric factors at ! C latitude DLAT. ! C DGR: resonance factor for gravimetric factors. ! C DHLAT: array(1..12) containing the Love-numbers h at ! C latitude DLAT. ! C DHR: resonance factor for the Love-number h(2,1). ! C DKLAT: array(1..12) containing the Love-numbers k at ! C latitude DLAT. ! C DKR: resonance factor for the Love-number k(2,1). ! C DLLAT: array(1..12) containing the Shida-numbers l at ! C latitude DLAT. ! C DLR: resonance factor for the Shida-number l(2,1). ! C DTLAT: array(1..12) containing the tilt factors at ! C latitude DLAT. ! C ! C Reference: ! C ---------- ! C ! C Dehant, V. (1987): Tidal parameters for an inelastic Earth. ! C Physics of the Earth and Planetary Interiors, 49, 97-116, ! C 1987. ! C ! C Wahr, J.M. (1981): Body tides on an elliptical, rotating, ! C elastic and oceanless earth. Geophysical Journal of the Royal ! C Astronomical Society, vol. 64, 677-703, 1981. ! C ! C Zschau, J. and R. Wang (1987): Imperfect elasticity in the ! C Earth's mantle. Implications for earth tides and long period ! C deformations. Proceedings of the 9th International Symposium ! C on Earth Tides, New York 1987, pp. 605-629, editor J.T. Kuo, ! C Schweizerbartsche Verlagsbuchhandlung, Stuttgart 1987. ! C ! C Routine creation: 1993.07.03 by Hans-Georg Wenzel, ! C Black Forest Observatory, ! C Universitaet Karlsruhe, ! C Englerstr. 7, ! C D-76128 KARLSRUHE, ! C Germany. ! C Tel: 0049-721-6082307, ! C FAX: 0049-721-694552. ! C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de ! C Last modification: 1996.05.25 by Hans-Georg Wenzel. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IMPLICIT DOUBLE PRECISION (D) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C The following DIMENSION statement is concerning the elastic ! C Earth model for the different degree and order constituents. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C COMMON /LOVE/ contains gravimeter factors, Love-numbers, Shida- ! C numbers and tilt factors for degree 2...4 at latitude DLAT: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C COMMON /CONST/: ! C DPI... 3.1415.... DPI2... 2.D0*DPI ! C DRAD... DPI/180.D0 DRO... 180.D0/DPI ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! COMMON /CONST/ DPI,DPI2,DRAD,DRO SAVE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C The following DATA statements are concerning the elastic ! C Earth model for the different degree and order constituents. ! C The latitude dependency is not given for all constituents in ! C the Wahr-Dehant-Zschau model !!!!!! ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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/ C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Shida-numbers: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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/ C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Definition of parameters of Geodetic Reference System 1980. ! C DEA is major semi axis in meter. ! C DEE is square of first excentricity (without dimnension). ! C DEGM is geocentric gravitational constant in m*3/s**2. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DATA DEA/6378137.00D0/,DEE/6.69438002290D-3/ C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Define resonance frequency and resonance factors: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DOMR=15.073729D0 DOM0=13.943036D0 DGR =-0.000625D0 DHR =-0.002505D0 DKR =-0.001261D0 DLR =0.0000781D0 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C DCLAT is cos and DSLAT is sin of ellipsoidal latitude. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DCLAT=DCOS(DLAT*DRAD) DSLAT=DSIN(DLAT*DRAD) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Compute ellipsoidal curvature radius DN in meter. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DN=DEA/DSQRT(1.D0-DEE*DSLAT**2) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Compute geocentric latitude DPSI in degree: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Compute latitude dependent gravimeter factors DG: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DO 110 I=1,12 110 DGLAT(I)=DG0(I)+DGP(I)*DLATP(I)+DGM(I)*DLATM(I) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Compute latitude dependent LOVE-numbers DH (for vertical ! C displacement): ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DO 120 I=1,12 120 DHLAT(I)=DH0(I)+DHP(I)*DLATP(I)+DHM(I)*DLATM(I) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Compute latitude dependent LOVE-numbers DK: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DO 130 I=1,12 130 DKLAT(I)=DK0(I)+DKP(I)*DLATP(I)+DKM(I)*DLATM(I) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Compute latitude dependent SHIDA-numbers DL: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DO 140 I=1,12 140 DLLAT(I)=DL0(I)+DLP(I)*DLATP(I)+DLM(I)*DLATM(I) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Compute latitude dependent tilt factors DT: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Print out of parameters: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! WRITE(IUN16,17001) DOM0,DOMR,DGR,DHR,DKR,DLR,DTR I=0 WRITE(IUN16,17002) DLAT DO 300 L=2,4 WRITE(IUN16,17004) DO 300 M=0,L I=I+1 WRITE(IUN16,17003) L,M,DGLAT(I),DHLAT(I),DKLAT(I),DLLAT(I), 1 DTLAT(I) 300 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Format statements: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 17001 FORMAT(/6x,'Routine ETLOVE, version 1996.05.25.'/ 1 6x,'Latitude dependent parameters for an elliptical, rotating,'/ 2 6x,'inelastic and oceanless Earth from Wahr-Dehant-Zschau model.' 3 // 4 6x,'frequency of wave O1:',F10.6,' deg per hour'/ 5 6x,'resonance frequency :',F10.6,' deg per hour'// 6 6x,'resonance factor for G:',F10.6/ 7 6x,'resonance factor for h:',F10.6/ 8 6x,'resonance factor for k:',F10.6/ 9 6x,'resonance factor for l:',F10.6/ * 6x,'resonance factor for T:',F10.6/) 17002 FORMAT(// 1 6x,'Latitude dependent elastic parameters'// 2 6x,'ellipsoidal latitude:',F10.4,' deg'// 3 6x,'G is gravimetric factor delta'/ 4 6x,'h is LOVE-number h'/ 5 6x,'k is LOVE-number k'/ 6 6x,'l is SHIDA-number l'/ 7 6x,'T is tilt factor gamma'// 8 6x,'degree order G h k l', 9' T') 17003 FORMAT(6x,2I7,5F10.6) 17004 FORMAT(' ') RETURN END C SUBROUTINE ETPHAS(IUN16,IPRINT,IMODEL,DLON,DJULD) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C ! C Routine ETPHAS, version 1996.08.03 Fortran 90. ! C ! C The routine ETPHAS computes phases and frequencies for the tidal ! C waves using different tidal potential catalogues which use ! C the Hartmann and Wenzel (1995) normalization. ! C ! C All variables with D as first character are DOUBLE PRECISION. ! C ! C Input parameter description: ! C ---------------------------- ! C ! C IUN16: Formatted line printer unit. ! C IPRINT: Printout parameter. ! C for IPRINT = 0, nothing will be printed. ! C for IPRINT = 1, a short list will be printed. ! C for IPRINT = 2, a long list will be printed ! C (including the tidal potential development). ! C IMODEL: Parameter for selecting the tidal potential ! C development. ! C IMODEL = 1: Doodson (1921) tidal potential develop- ! C ment with 378 waves. ! C IMODEL = 2: Cartwright-Taylor-Edden (1973) tidal ! C potential development with 505 waves. ! C IMODEL = 3: Buellesfeld (1985) tidal potential ! C development with 656 waves. ! C IMODEL = 4: Tamura (1987) tidal potential develop- ! C ment with 1200 waves. ! C IMODEL = 5: Xi (1989) tidal potential catalogue ! C 2933 waves. ! C IMODEL = 6: Roosbeek (1995) tidal potential ! C catalogue with ?? waves. ! C IMODEL = 7: Hartmann and Wenzel (1995) tidal ! C potential catalogue with 12935 waves. ! C DLON: Ellipsoidal longitude referring to Geodetic ! C Reference System 1980 in degree, positive east of ! C Greenwhich. ! C DJULD: Julian date of the initial epoch of tidal force ! C development. ! C ! C Output parameter description: ! C ----------------------------- ! C ! C There are no output parameters. The computes phases are trans- ! C to the calling program unit by COMMON /TIDWAVE/. ! C ! C COMMON /TIDWAVE/: contains tidal waves ! C ! C NW: Number of defined tidal waves. ! C IWNR: INTEGER array (1:12935) of wave numbers. ! C IAARG: INTEGER array (1:12935,1:12) of astronomical ! C argument numbers. ! C DX0: DOUBLE PRECISION array (1:12935) of cos-coeffi- ! C cients of the tidal component in units of the tidal ! C component. ! C DX1: DOUBLE PRECISION array (1:12935) of time deriva- ! C tives of cos-coefficients of the tidal component. ! C DY0: DOUBLE PRECISION array (1:12935) of sin-coeffi- ! C cients of the tidal component in units of the tidal ! C component. ! C DY1: DOUBLE PRECISION array (1:12935) of time deriva- ! C tives of sin-coefficients of the tidal component. ! C ! C component unit of unit of ! C IC DX0,DY0 DX1,DY1 ! C -1 m**2/s**2 m**2/s**2 per Julian century ! C 0 nm/s**2 nm/s**2 per Julina century ! C 1 mas mas per Julian century ! C 2 mm mm per Julian century ! C 3 mm mm per Julian century ! C 4 nstr nstr per Julian cenrury ! C 5 nstr nstr per Julian century ! C 6 nstr nstr per Julian century ! C 7 nstr nstr per Julian century ! C 8 nstr nstr per Julian century ! C 9 mm mm per Julian century ! C ! C DTHPH: DOUBLE PRECISION array (1:12935) of tidal phases ! C in radians at initial epoch. ! C DTHFR: DOUBLE PRECISION array (1:12935) of tidal ! C frequencies in radian per hour. ! C DBODY: DOUBLE PRECISION array (1:12935) of body tide ! C amplitude factors for tidal gravity and tidal tilt. ! C In order to compute the body tide, the coefficients ! C DX0, DX1, DY0 and DY1 have to be multiplied by ! C DBODY. ! C ! C Used routines: ! C -------------- ! C ! C ETASTN: computes astronomical elements. ! C ETJULN: computes Julian date. ! C ETDDTA: computes the difference TDT minus UTC (called by ETASTN).! C ETPOLC: computes the difference DUT1 = UT1 - UTC. ! C ! C Numerical accuracy: ! C ------------------- ! C ! C The routine has been tested under operation systems UNIX and ! C MS-DOS with 15 digits in DOUBLE PRECISION. ! C ! C Routine creation: 1988.04.27 by Hans-Georg Wenzel, ! C Black Forest Observatory, ! C Universitaet Karlsruhe, ! C Englerstr. 7, ! C D-76128 KARLSRUHE 1, ! C Germany. ! C Tel: 0049-721-6082307, ! C FAX: 0049-721-694552. ! C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de ! C Last modification: 1996.08.04 by Hans-Georg Wenzel. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IMPLICIT DOUBLE PRECISION (D) IMPLICIT INTEGER (I-N) CHARACTER CMODEL(7)*20 DOUBLE PRECISION DAS(11),DASP(11) COMMON /TIDPHAS/ DPK(25) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C The following DIMENSION statement is concerning the number of ! C waves of the tidal potential development, which is 12935. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! COMMON /TIDWAVE/ NW,IWNR(12935),IAARG(12935,12),DX0(12935), 1 DX1(12935),DY0(12935),DY1(12935),DTHPH(12935),DTHFR(12935), 2 DBODY(12935) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C COMMON /CONST/: ! C DPI... 3.1415.... ! C DPI2... 2.D0*DPI ! C DRAD... DPI/180.D0 ! C DRO... 180.D0/DPI ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! COMMON /CONST/ DPI,DPI2,DRAD,DRO SAVE DATA IUN30/30/,IUN31/31/ DATA CMODEL/'Doodson 1921 ', 1 'CTED 1973 ','Buellesfeld 1985 ', 2 'Tamura 1987 ','Xi 1989 ', 3 'Roosbeek 1995 ','Hartmann+Wenzel 1995'/ IF(IPRINT.GT.0) WRITE(IUN16,17001) CMODEL(IMODEL) 1000 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Interpolate DUT1: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CALL ETPOLC(IUN16,IUN30,IUN31,IPRINT,DJULD,DCLAT,DSLAT, 1 DCLON,DSLON,DPOLX,DPOLY,DUT1,DTAI,DLOD,DGPOL,DGPOLP,DGLOD,NERR) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Compute astronomical elements for initial epoch: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CALL ETASTN(IUN16,IPRINT,IMODEL,DLON,DJULD,DUT1,DAS,DASP,DDT0) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Compute phases and frequencies: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DO 1110 IW=1,NW DC2=0.D0 DC3=0.D0 DO 1140 J=1,11 DC2=DC2+DBLE(IAARG(IW,J))*DAS(J) 1140 DC3=DC3+DBLE(IAARG(IW,J))*DASP(J) LI=IAARG(IW,12) JCOF=(LI+1)*LI/2-2+IAARG(IW,1) DC2=DC2+DPK(JCOF) 1160 DC2=DMOD(DC2,360.D0) IF(DC2.GE.0.D0) GOTO 1170 DC2=DC2+360.D0 GOTO 1160 1170 DTHPH(IW)=DC2*DRAD DTHFR(IW)=DC3*DRAD 1110 CONTINUE IF(IPRINT.EQ.0) RETURN WRITE(IUN16,17002) NW WRITE(IUN16,17003) RETURN C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Format statements: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 17001 FORMAT(' Routine ETPHAS, version 1996.08.04.'// 1' Tidal component development from tidal potential development.'// 2 1X,A13,' tidal potential development is used.'/) 17002 FORMAT(//' Routine ETPHAS, version 1996.08.04.'/ 1'New phases and frequencies computes for',I6,' waves.') 17003 FORMAT(///' ***** Routine ETPHAS finished execution.'/) END C SUBROUTINE ETPOLC(IUN16,IUN30,IUN31,IPRINT,DJULD,DCLAT,DSLAT, 1 DCLON,DSLON,DPOLX,DPOLY,DUT1,DTAI,DLOD,DGPOL,DGPOLP,DGLOD,NERR) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C ! C Routine ETPOLC, version 1996.05.25 Fortran 90. ! C ! C The routine ETPOLC returns pole coordinates and correction DUT1 ! C read from either formatted file on IUN30 or unformatted direct ! C access file on IUN31. In case that direct access file IUN31 does ! C not exist, it will be established by routine ETPOLC with file ! C /home/hwz/eterna34/commdat/etpolut1.uft. ! C ! C All variables with D as first character are DOUBLE PRECISION. ! C ! C Input parameter description: ! C ---------------------------- ! C ! C IUN16: Formatted line printer unit. ! C IUN30: Formatted file containing pole coordinates, DUT1 ! C and DTAI (e.g. file etpolut1.dat). ! C IUN31: Unformatted direct access file containing pole ! C coordinates, DUT1 and DTAI. This file will be ! C opened as file etpolut1.uft during the execution of ! C routine ETPOLC with STATUS=OLD if it exists and ! C with STATUS=NEW, if it does not exist. If the ! C file does not yet exist, etpolut1.uft will be ! C established during the execution of routine ! C ETPOLC. ! C IPRINT: Printout parameter. ! C for IPRINT = 0, nothing will be printed. ! C for IPRINT = 1, a short list will be printed. ! C for IPRINT = 2, a long list will be printed ! C DJULD: Julian date of the epoch, for which pole ! C coordinates, DUT1 and DTAI will be returned. ! C DCLAT: COS of latitude. ! C DSLAT: SIN of latitude. ! C DCLON: COS of longitude. ! C DSLON: SIN of longitude. ! C ! C Output parameter description: ! C ----------------------------- ! C ! C DPOLX: X-pole coordinate in arc sec. ! C DPOLY: Y-pole coordinate in arc sec. ! C DUT1: Difference UT1 minus UTC in sec. ! C DTAI: Difference TAI minus UT1 in sec. ! C DLOD: Length of day - 86400 sec in sec. ! C DGPOL: Pole tide in nm/s**2 for a rigid earth. ! C DGPOLP: Time derivative of pole tide for a rigid earth in ! C nm/s**2 per day. ! C DGLOD: Gravity variation due to variation of the earth's ! C rotation in nm/s**2. ! C NERR: Error code, counts the number of errors which ! C happened during the actual call of routine ETPOLC. ! C For NERR > 0, the output parameters DPOLX, DPOLY, ! C DUT1, DTAI, DLOD, DGPOL, DGPOLP do not contain ! C valid information (all set to zero). ! C For NERR=0, the output parameters DPOLX, DPOLY, ! C DUT1 and DTAI contain valid information. ! C ! C Execution time: ! C --------------- ! C ! C 3.02 microsec per call on a 100 MHz Pentium using Lahey LF90 ! C compiler if the file ETPOLC.UFT exists already. ! C ! C Routine creation: 1993.08.31 by Hans-Georg Wenzel, ! C Black Forest Observatory, ! C Universitaet Karlsruhe, ! C Englerstr. 7, ! C D-76128 KARLSRUHE, ! C Germany. ! C Tel: 0049-721-6082307, ! C FAX: 0049-721-694552. ! C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de ! C Last modification: 1996.05.25 by Hans-Georg Wenzel. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IMPLICIT DOUBLE PRECISION (D) IMPLICIT INTEGER (I-N) CHARACTER CHEAD(8)*10,CENDH*10 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C COMMON /CONST/: to be initialzed by BLOCK DATA before the call ! C of routine ETPOLC: ! C DPI... 3.1415.... ! C DPI2... 2.D0*DPI ! C DRAD... DPI/180.D0 ! C DRO... 180.D0/DPI ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! COMMON /CONST/ DPI,DPI2,DRAD,DRO SAVE DATA DOM/7.292115D-5/,DA/6378137.D0/ DATA CENDH/'C*********'/,ISTART/1/,IMJDO/0/ NERR=0 IF(ISTART.EQ.0) GOTO 1000 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Test, whether there exist already unformatted file ETPOLUT1.UFT: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C HP-UNIX: OPEN(UNIT=IUN31,FILE='../commdat/etpolut1.uft', C HP-UNIX: 1 FORM='UNFORMATTED',STATUS='OLD',ACCESS='DIRECT',RECL=80,ERR=11) C MS-DOS: OPEN(UNIT=IUN31,FILE='/home/hwz/eterna34/commdat/etpolut2.uft', 1 FORM='UNFORMATTED',STATUS='OLD',ACCESS='DIRECT',RECL=32,ERR=11) WRITE(*,'(A$)')' FILE etpolut2.uft is existing ' READ(IUN31,REC=1) IFIRST,ILAST write(*,*)ifirst,ilast ISTART=0 GOTO 1000 C HP-UNIX: 11 OPEN(UNIT=IUN31,FILE='../commdat/etpolut2.uft', C HP-UNIX: 1 FORM='UNFORMATTED',STATUS='NEW',ACCESS='DIRECT',RECL=80) C MS-DOS: 11 OPEN(UNIT=IUN31,FILE='/home/hwz/eterna34/commdat/etpolut2.uft', 1 FORM='UNFORMATTED',STATUS='NEW',ACCESS='DIRECT',RECL=32) c REWIND IUN31 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Read file header of tidal potential file on unit IUN30: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(IPRINT.EQ.0) GOTO 10 WRITE(IUN16,17001) 10 CONTINUE READ(IUN30,17002) (CHEAD(I),I=1,8) WRITE(IUN16,17003) (CHEAD(I),I=1,8) 100 READ(IUN30,17002) (CHEAD(I),I=1,8) IF(IPRINT.GT.1) WRITE(IUN16,17003) (CHEAD(I),I=1,8) IF(CHEAD(1).NE.CENDH) GOTO 100 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Read data: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IREC=2 ILAST=0 200 READ(IUN30,17004) IDAT,ITIM,DMODJI,DPOLX,DPOLY,DUT1,DTAI IF(IDAT.EQ.99999999) GOTO 300 IF(IREC.EQ.2) IFIRST=DMODJI WRITE(IUN31,REC=IREC) DPOLX,DPOLY,DUT1,DTAI IF(IPRINT.GT.1) WRITE(IUN16,17005) IDAT,ITIM,IREC,DMODJI,DPOLX, 1 DPOLY,DUT1,DTAI ILAST=IREC IREC=IREC+1 GOTO 200 300 CONTINUE WRITE(IUN31,REC=1) IFIRST,ILAST write(*,*)ifirst,ilast ISTART=0 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Read pole coordinates, DUT1 and DTAI from direct access unit ! C IUN31: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1000 DMODJD=DJULD-2400000.5D0 IMJD=DMODJD C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C DT is time difference referring to central sample point in days: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DT=DMODJD-DBLE(IMJD) DT2=DT*DT IREC=IMJD-IFIRST+2 IF(IREC.LT.2) THEN DPOLX=0.D0 DPOLY=0.D0 DUT1 =0.D0 DTAI =0.D0 DLOD =0.D0 DGPOL=0.D0 DGPOLP=0.D0 NERR=1 RETURN ENDIF IF(IREC.GT.ILAST-1) THEN C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Use pole coordinates and DUT1 from last tabulated day: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! READ(IUN31,REC=ILAST) DPOLX,DPOLY,DUT1,DTAI DLOD =0.D0 DGPOL=DOM**2*DA*2.D0*DCLAT*DSLAT*(DPOLX*DCLON-DPOLY*DSLON)* 1 DRAD/3600.D0*1.D9 DGPOLP=0.D0 NERR=1 RETURN ENDIF IF(IMJD.EQ.IMJDO) GOTO 1100 READ(IUN31,REC=IREC-1) DPOLX1,DPOLY1,DUT12,DTAI1 READ(IUN31,REC=IREC) DPOLX2,DPOLY2,DUT12,DTAI2 READ(IUN31,REC=IREC+1) DPOLX3,DPOLY3,DUT13,DTAI3 IMJDO=IMJD C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Quadratic interpolation for pole coordinates and DTAI: ! C Linear interpolation for DUT1: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1100 DPOLXA0=DPOLX2 DPOLXA1=(DPOLX3-DPOLX1)*0.5D0 DPOLXA2=(DPOLX1-2.D0*DPOLX2+DPOLX3)*0.5D0 C DPOLYA0=DPOLY2 DPOLYA1=(DPOLY3-DPOLY1)*0.5D0 DPOLYA2=(DPOLY1-2.D0*DPOLY2+DPOLY3)*0.5D0 C DTAIA0=DTAI2 DTAIA1=(DTAI3-DTAI1)*0.5D0 DTAIA2=(DTAI1-2.D0*DTAI2+DTAI3)*0.5D0 C DUT10=DUT12 DDUT1=DUT13-DUT12 IF(DDUT1.GT. 0.9D0) DDUT1=DDUT1-1.D0 IF(DDUT1.LT.-0.9D0) DDUT1=DDUT1+1.D0 DLOD = DTAIA1+2.D0*DTAIA2*DT DGLOD=2.D0*DLOD*DOM**2*DA*DCLAT*DCLAT*1.D9/86400.D0 C DPOLX=DPOLXA0+DT*DPOLXA1+DT2*DPOLXA2 DPOLY=DPOLYA0+DT*DPOLYA1+DT2*DPOLYA2 DUT1 =DUT10 +DT*DDUT1 DTAI =DTAIA0 +DT*DTAIA1 +DT2*DTAIA2 C DGPOL=DOM**2*DA*2.D0*DCLAT*DSLAT*(DPOLX*DCLON-DPOLY*DSLON)* 1 DRAD/3600.D0*1.D9 DPOLXP=DPOLXA1+2.D0*DPOLXA2*DT DPOLYP=DPOLYA1+2.D0*DPOLYA2*DT DGPOLP=DOM**2*DA*2.D0*DCLAT*DSLAT*(DPOLXP*DCLON-DPOLYP*DSLON)* 1 DRAD/3600.D0*1.D9 RETURN C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Format statements: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 17001 FORMAT(' Routine ETPOLC, version 1996.05.25.'// 1' Pole coordinates, DUT1, DTAI and pole tides from IERS data.'//) 17002 FORMAT(8A10) 17003 FORMAT(1X,8A10) 17004 FORMAT(I8,1X,I6,F10.3,5F10.5) 17005 FORMAT(I9,1X,2I6,F10.3,5F10.5) END C SUBROUTINE ETPOTS(IUN14,IUN16,IUN24,IPRINT,IMODEL,DLAT,DLON,DH, 1 DGRAV,DAZ,IC,DJULD,DAMIN) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C ! C Routine ETPOTS, version 1996.08.05 Fortran 90. ! C ! C The routine ETPOTS computes amplitudes, phases, frequencies and ! C body tide amplitude factors for a number of different Earth tide ! C components using different tidal potential catalogues which use ! C the Hartmann and Wenzel (1995) normalization. ! C ! C Attention: This routine has finally not been tested for vertical ! C and horizontal displacements and for shear tidal ! C strain !!!! ! C ! C All variables with D as first character are DOUBLE PRECISION. ! C ! C Input parameter description: ! C ---------------------------- ! C ! C IUN14: Formatted unit, on which the tidal potential ! C development has to be stored before the execution ! C of routine ETPOTS (e.g. file hw95s.dat). ! C IUN16: Formatted line printer unit. ! C IUN24: Unformatted copy of IUN14. This unit will be opened ! C e.g. as file hw95s.uft during the execution of ! C routine ETPOTS with STATUS=OLD if it exists and ! C with STATUS=NEW, if it does not exist. If the file ! C does not yet exist, it will be established during ! C the execution of routine ETPOTS. ! C IPRINT: Printout parameter. ! C for IPRINT = 0, nothing will be printed. ! C for IPRINT = 1, a short list will be printed. ! C for IPRINT = 2, a long list will be printed ! C (including the tidal potential development). ! C IMODEL: Parameter for selecting the tidal potential ! C development. ! C IMODEL = 1: Doodson (1921) tidal potential develop- ! C ment with 378 waves. ! C IMODEL = 2: Cartwright-Taylor-Edden (1973) tidal ! C potential development with 505 waves. ! C IMODEL = 3: Buellesfeld (1985) tidal potential ! C development with 656 waves. ! C IMODEL = 4: Tamura (1987) tidal potential develop- ! C ment with 1200 waves. ! C IMODEL = 5: Xi (1989) tidal potential catalogue ! C 2933 waves. ! C IMODEL = 6: Roosbeek (1995) tidal potential ! C catalogue with ?? waves. ! C IMODEL = 7: Hartmann and Wenzel (1995) tidal ! C potential catalogue with 12935 waves. ! C DLAT: Ellipsoidal latitude referring to Geodetic ! C Reference System 1980 in degree. ! C DLON: Ellipsoidal longitude referring to Geodetic ! C Reference System 1980 in degree, positive east of ! C Greenwhich. ! C DH: Ellipsoidal height referring to Geodetic Reference ! C System 1980 in meter. ! C DGRAV: Gravity in m/s**2. If the gravity is input below ! C 1 m/s**2, the gravity will be replaced by the ! C computed normal gravity for reference system GRS80. ! C DAZ: Azimuth in degree from north direction (only valid ! C for tidal tilt, horizontal displacement, and ! C horizontal strain). ! C IC: Earth tide component to be computed. ! C IC=-1: tidal potential in m**2/s**2. ! C IC= 0: vertical tidal acceleration (gravity tide), ! C in nm/s**2 (positive downwards). ! C IC= 1: horizontal tidal acceleration (tidal tilt) ! C in azimuth DAZ in mas = arc sec/1000. ! C IC= 2: vertical tidal displacement, geodetic ! C coefficients in mm (positive upwards). ! C IC= 3: horizontal tidal displacement in azimuth ! C DAZ in mm. ! C IC= 4: vertical tidal strain in 10**-9 = nstr. ! C IC= 5: horizontal tidal strain in azimuth DAZ ! C in 10**-9 = nstr. ! C IC= 6: areal tidal strain in 10**-9 = nstr. ! C IC= 7: shear tidal strain in 10**-9 = nstr. ! C IC= 8: volume tidal strain in 10**-9 = nstr. ! C IC= 9: ocean tides, geodetic coefficients in ! C millimeter. ! C DJULD: Julian date of the initial epoch of tidal force ! C development. ! C DAMIN: Truncation parameter for the amplitude of tidal ! C waves to be used in m**2/s**2. Only tidal waves ! C with amplitudes greater or equal DAMIN will be ! C used. ! C ! C Rms error of gravity tides compited from HW95 tidal ! C potential catalogue versus amaplitude threshold, ! C as computed from comparison with benchmark gravity ! C tide series BFDE403A ! C ! C DAMIN no. of rms error min. error max.error ! C [m**2/s**2] waves [nm/s**2] [nm/s**2] [nm/s**2] ! C ! C 1.00*10**-1 11 88.403330 -321.492678 297.866988 ! C 3.16*10**-2 28 27.319455 -108.174675 109.525103 ! C 1.00*10**-2 45 14.449139 -62.286861 67.322802 ! C 3.16*10**-3 85 6.020159 -32.560229 28.931931 ! C 1.00*10**-3 158 2.249690 -14.587415 11.931120 ! C 3.16*10**-4 268 0.978419 -6.780051 5.934767 ! C 1.00*10**-4 441 0.436992 -3.049676 2.943019 ! C 3.16*10**-5 768 0.173071 -1.331572 1.242490 ! C 1.00*10**-5 1 273 0.068262 -0.520909 0.484510 ! C 3.16*10**-6 2 052 0.029229 -0.217114 0.229504 ! C 1.00*10**-6 3 359 0.011528 -0.099736 0.085920 ! C 3.16*10**-7 5 363 0.004706 -0.038247 0.035942 ! C 1.00*10**-7 8 074 0.001999 -0.019407 0.017684 ! C 3.16*10**-8 10 670 0.001391 -0.012350 0.012287 ! C 1.00*10**-8 12 234 0.001321 -0.010875 0.011307 ! C ! C Output parameter description: ! C ----------------------------- ! C ! C There are no output parameters. The computed arrays are trans- ! C ferred to the calling program unit by COMMON /TIDWAVE/. ! C ! C COMMON /TIDWAVE/: contains tidal waves ! C ! C NW: Number of defined tidal waves. ! C IWNR: INTEGER array (1:12935) of wave numbers. ! C IAARG: INTEGER array (1:12935,1:12) of astronomical ! C argument numbers. ! C DX0: DOUBLE PRECISION array (1:12935) of cos-coeffi- ! C cients of the tidal component in units of the tidal ! C component. ! C DX1: DOUBLE PRECISION array (1:12935) of time deriva- ! C tives of cos-coefficients of the tidal component. ! C DY0: DOUBLE PRECISION array (1:12935) of sin-coeffi- ! C cients of the tidal component in units of the tidal ! C component. ! C DY1: DOUBLE PRECISION array (1:12935) of time deriva- ! C tives of sin-coefficients of the tidal component. ! C ! C component unit of unit of ! C IC DX0,DY0 DX1,DY1 ! C -1 m**2/s**2 m**2/s**2 per Julian century ! C 0 nm/s**2 nm/s**2 per Julina century ! C 1 mas mas per Julian century ! C 2 mm mm per Julian century ! C 3 mm mm per Julian century ! C 4 nstr nstr per Julian cenrury ! C 5 nstr nstr per Julian century ! C 6 nstr nstr per Julian century ! C 7 nstr nstr per Julian century ! C 8 nstr nstr per Julian century ! C 9 mm mm per Julian century ! C ! C DTHPH: DOUBLE PRECISION array (1:12935) of tidal phases ! C in radians at initial epoch. ! C DTHFR: DOUBLE PRECISION array (1:12935) of tidal ! C frequencies in radian per hour. ! C DBODY: DOUBLE PRECISION array (1:12935) of body tide ! C amplitude factors for tidal gravity and tidal tilt. ! C In order to compute the body tide, the coefficients ! C DX0, DX1, DY0 and DY1 have to be multiplied by ! C DBODY. ! C ! C Used routines: ! C -------------- ! C ! C ETASTN: computes astronomical elements. ! C ETGCON: computes geodetic coefficients. ! C ETJULN: computes Julian date. ! C ETLOVE: computes latitude dependent elastic parameters (called ! C ETGCOF). ! C ETDDTA: computes the difference TDT minus UTC (called by ETASTN).! C ETPOLC: computes the difference DUT1 = UT1 - UTC. ! C ! C Numerical accuracy: ! C ------------------- ! C ! C The routine has been tested under operation systems UNIX and ! C MS-DOS with 15 digits in DOUBLE PRECISION. ! C ! C Routine creation: 1988.04.27 by Hans-Georg Wenzel, ! C Black Forest Observatory, ! C Universitaet Karlsruhe, ! C Englerstr. 7, ! C D-76128 KARLSRUHE 1, ! C Germany. ! C Tel: 0049-721-6082307, ! C FAX: 0049-721-694552. ! C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de ! C Last modification: 1996.08.05 by Hans-Georg Wenzel. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IMPLICIT DOUBLE PRECISION (D) IMPLICIT INTEGER (I-N) LOGICAL LEX24 CHARACTER CHEAD(8)*10,CENDH*10,CUNIT(11)*8 CHARACTER CMODEL(7)*20,CFFILE(7)*64,CUFILE(7)*64 CHARACTER CBOD*2,CWAVE*4 INTEGER NS(11) DOUBLE PRECISION DAS(11),DASP(11),DGK(25),DPK(25) COMMON /TIDPHAS/ DPK C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C The following DIMENSION statement is concerning the number of ! C waves of the tidal potential development, which is 12935. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! COMMON /TIDWAVE/ NW,IWNR(12935),IAARG(12935,12),DX0(12935), 1 DX1(12935),DY0(12935),DY1(12935),DTHPH(12935),DTHFR(12935), 2 DBODY(12935) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C The following DIMENSION statement is concerning the elastic ! C Earth model for the different degree and order constituents. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DOUBLE PRECISION DELTA(25) COMMON /UNITS/ CUNIT,IC2 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C COMMON /CONST/: ! C DPI... 3.1415.... ! C DPI2... 2.D0*DPI ! C DRAD... DPI/180.D0 ! C DRO... 180.D0/DPI ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! COMMON /CONST/ DPI,DPI2,DRAD,DRO C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C COMMON /LOVE/ contains gravimeter factors, LOVE-numbers, SHIDA- ! C numbers and tilt factors for degree 2...4 at latitude DLAT: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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 SAVE DATA MAXNW/12935/ DATA CENDH/'C*********'/ DATA IUN30/30/,IUN31/31/ DATA CMODEL/'Doodson 1921 ', 1 'CTED 1973 ','Buellesfeld 1985 ', 2 'Tamura 1987 ','Xi 1989 ', 3 'Roosbeek 1995 ','Hartmann+Wenzel 1995'/ C HP-UNIX: C HP-UNIX: DATA CFFILE/ '../commdat/doodsehw.dat', C HP-UNIX: 2 '../commdat/cted73hw.dat', C HP-UNIX: 3 '../commdat/buellehw.dat', C HP-UNIX: 4 '../commdat/tamurahw.dat', C HP-UNIX: 5 '../commdat/xi1989hw.dat', C HP-UNIX: 6 '../commdat/ratgp95.dat', C HP-UNIX: 7 '../commdat/hw95s.dat'/ C HP-UNIX: DATA CUFILE/ '../commdat/doodsehw.uft', C HP-UNIX: 2 '../commdat/cted73hw.uft', C HP-UNIX: 3 '../commdat/buellehw.uft', C HP-UNIX: 4 '../commdat/tamurahw.uft', C HP-UNIX: 5 '../commdat/xi1989hw.uft', C HP-UNIX: 6 '../commdat/ratgp95.uft', C HP-UNIX: 7 '../commdat/hw95s.uft'/ C MS-DOS: DATA CFFILE/ '/home/hwz/eterna34/commdat/doodsehw.dat', 2 '/home/hwz/eterna34/commdat/cted73hw.dat', 3 '/home/hwz/eterna34/commdat/buellehw.dat', 4 '/home/hwz/eterna34/commdat/tamurahw.dat', 5 '/home/hwz/eterna34/commdat/xi1989hw.dat', 6 '/home/hwz/eterna34/commdat/ratgp95.dat', 7 '/home/hwz/eterna34/commdat/hw95s.dat'/ DATA CUFILE/ '/home/hwz/eterna34/commdat/doodseh2.uft', 2 '/home/hwz/eterna34/commdat/cted73h2.uft', 3 '/home/hwz/eterna34/commdat/buelleh2.uft', 4 '/home/hwz/eterna34/commdat/tamurah2.uft', 5 '/home/hwz/eterna34/commdat/xi1989h2.uft', 6 '/home/hwz/eterna34/commdat/ratgp952.uft', 7 '/home/hwz/eterna34/commdat/hw95s2.uft'/ IF(IPRINT.GT.0) WRITE(IUN16,17001) CMODEL(IMODEL) OPEN(UNIT=IUN14,FILE=CFFILE(IMODEL),FORM='FORMATTED', 1 STATUS='OLD') REWIND(IUN14) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Test, whether there exist already the unformatted tidal ! C potential catalogue file: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! OPEN(UNIT=IUN24,FILE=CUFILE(IMODEL),FORM='UNFORMATTED', 1 STATUS='OLD',ERR=11) LEX24=.TRUE. WRITE(*,17002)CUFILE(IMODEL) REWIND IUN24 GOTO 12 11 OPEN(UNIT=IUN24,FILE=CUFILE(IMODEL),FORM='UNFORMATTED', 1 STATUS='NEW') LEX24=.FALSE. REWIND IUN14 12 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Compute geodetic coefficients and body tide amplitude factors ! C for the WAHR-DEHANT-ZSCHAU model. The NDFW resonance is ! C approximated by ! C ! C G0 - GR*(DOM - DOM0)/(DOMR - DOM), ! C ! C similar equations hold for the other components. ! C ! C Gravimetric amplitude factors, LOVE numbers h and k for zero to ! C third degree tidal potential have been taken from DEHANT 1987, ! C table 7, 8 and 9 for elliptical, uniformly rotating, oceanless ! C Earth with liquid outer core and inelastic mantle (PREM Earth ! C model with inelastic mantle from ZSCHAU) and for the fourth ! C degree from DEHANT et al. 1989, table 6). The resonance factors ! C GR have been computed to fit the difference between body tide ! C amplitude factors at waves O1 and PSI1 from DEHANT 1987, PREM ! C model with elastic mantle (table 1...3). The NDFW resonance ! C frequency is 15.073729 degree per hour = 1.004915267 CPD UT, ! C taken from WAHR 1981 (because it is not given in any of DEHANT's ! C papers). ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CALL ETGCON(IUN16,IPRINT,DLAT,DLON,DH,DGRAV,DAZ,IC,DGK,DPK) IC2=IC+2 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Define default body tide amplitude factors for components ! C IC=2...9. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DO 50 I=1,25 50 DELTA(I)=1.D0 DELTAR=0.D0 GOTO (100,200,300),IC2 GOTO 1000 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C IC=-1, compute body tide amplitude factors for tidal potential: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 100 CONTINUE DO 110 I=1,12 110 DELTA(I)=DKLAT(I) DELTAR=DKR GOTO 1000 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C IC=0, compute body tide amplitude factors for vertical component ! C (gravity tides): ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 200 CONTINUE DO 210 I=1,12 210 DELTA(I)=DGLAT(I) DELTAR=DGR GOTO 1000 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C IC=1: compute body tide amplitude factors for horizontal ! C component (tidal tilt): ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 300 CONTINUE DO 310 I=1,12 310 DELTA(I)=DTLAT(I) DELTAR=DKR-DHR 1000 CONTINUE DT2000=(DJULD-2451544.D0)/36525.0D0 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Interpolate DUT1: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CALL ETPOLC(IUN16,IUN30,IUN31,IPRINT,DJULD,DCLAT,DSLAT, 1 DCLON,DSLON,DPOLX,DPOLY,DUT1,DTAI,DLOD,DGPOL,DGPOLP,DGLOD,NERR) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Compute astronomical elements for initial epoch: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CALL ETASTN(IUN16,IPRINT,IMODEL,DLON,DJULD,DUT1,DAS,DASP,DDT0) IC2=IC+2 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Read file header of tidal potential file on unit IUN14: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(LEX24) THEN READ(IUN24) (CHEAD(I),I=1,8) ELSE READ(IUN14,17028) (CHEAD(I),I=1,8) WRITE(IUN24) (CHEAD(I),I=1,8) ENDIF WRITE(IUN16,17029) (CHEAD(I),I=1,8) 1100 CONTINUE IF(LEX24) THEN READ(IUN24) (CHEAD(I),I=1,8) ELSE READ(IUN14,17028) (CHEAD(I),I=1,8) WRITE(IUN24) (CHEAD(I),I=1,8) ENDIF IF(IPRINT.EQ.2) WRITE(IUN16,17029) (CHEAD(I),I=1,8) IF(CHEAD(1).NE.CENDH) GOTO 1100 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Compute tidal development for the specific component from tidal ! C potential development: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IW=1 NWFILE=0 NAMPL=0 NTRUNC=0 1110 CONTINUE 1120 CONTINUE C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Read tidal potential catalogue either from formatted or from ! C unformatted file. The format of the files is described in ! C Hartmann and Wenzel (1995a). ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(LEX24) THEN READ(IUN24) NRI,CBOD,LI,(NS(J),J=1,11),DFR,DC0I,DS0I,DC1I, 1 DS1I,CWAVE ELSE READ(IUN14,17006,END=2000) NRI,CBOD,LI,(NS(J),J=1,11),DFR, 1 DC0I,DS0I,DC1I,DS1I,CWAVE WRITE(IUN24) NRI,CBOD,LI,(NS(J),J=1,11),DFR,DC0I,DS0I,DC1I, 1 DS1I,CWAVE ENDIF IF(NRI.GT.MAXNW) GOTO 2000 NWFILE=NWFILE+1 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Truncation of the tidal potential catalogue: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DAM=DSQRT(DC0I**2+DS0I**2)*1.D-10 IF(DAM.LT.DAMIN) THEN NTRUNC=NTRUNC+1 GOTO 1110 ENDIF IF(IW.EQ.1) GOTO 1130 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Check if the astronomical arguments are identical to those of ! C the last stored wave (for Hartmann and Wenzel 1995 potential): ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IDIFF=(LI-IAARG(IW-1,12))**2 DO 1125 J=1,11 1125 IDIFF=IDIFF+(NS(J)-IAARG(IW-1,J))**2 IF(IDIFF.GT.0) GOTO 1130 C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Astronomical arguments are identical to those of last stored ! C wave. We will add up the coefficients for these two waves: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(IW-1.GT.1) IWNR(IW-1)=NRI JCOF=(LI+1)*LI/2-2+NS(1) DX0(IW-1)=DX0(IW-1)+DC0I*DGK(JCOF)*1.D-10 DY0(IW-1)=DY0(IW-1)+DS0I*DGK(JCOF)*1.D-10 DX1(IW-1)=DX1(IW-1)+DC1I*DGK(JCOF)*1.D-10 DY1(IW-1)=DY1(IW-1)+DS1I*DGK(JCOF)*1.D-10 GOTO 1110 1130 NAMPL=NAMPL+1 DC2=0.D0 DC3=0.D0 IAARG(IW,12)=LI DO 1140 J=1,11 IAARG(IW,J)=NS(J) DC2=DC2+DBLE(NS(J))*DAS(J) 1140 DC3=DC3+DBLE(NS(J))*DASP(J) JCOF=(LI+1)*LI/2-2+NS(1) DC2=DC2+DPK(JCOF) IWNR(IW)=NRI DX0(IW)=DC0I*DGK(JCOF)*1.D-10 DY0(IW)=DS0I*DGK(JCOF)*1.D-10 DX1(IW)=DC1I*DGK(JCOF)*1.D-10 DY1(IW)=DS1I*DGK(JCOF)*1.D-10 DBODY(IW)=DELTA(JCOF) IF(JCOF.EQ.2) DBODY(IW)=DELTA(JCOF)+DELTAR*(DC3-DOM0)/(DOMR-DC3) 1160 DC2=DMOD(DC2,360.D0) IF(DC2.GE.0.D0) GOTO 1170 DC2=DC2+360.D0 GOTO 1160 1170 CONTINUE DTHPH(IW)=DC2*DRAD DTHFR(IW)=DC3*DRAD IF(IPRINT.EQ.2) THEN DXTI=DX0(IW)+DX1(IW)*DT2000 DYTI=DY0(IW)+DY1(IW)*DT2000 DTHAM=DSQRT(DXTI**2+DYTI**2) WRITE(IUN16,17011) IW,CBOD,LI,NS(1),DTHAM,DC2,DC3,CWAVE, 1 DBODY(IW) ENDIF IW=IW+1 IF(IW.GT.MAXNW) GOTO 5000 GOTO 1110 2000 CONTINUE NW=IW-1 CLOSE(IUN14) IF(IPRINT.EQ.0) RETURN WRITE(IUN16,17010) NWFILE,NTRUNC,NW WRITE(IUN16,17030) RETURN 5000 CONTINUE WRITE(IUN16,17050) NW,MAXNW STOP C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Format statements: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 17001 FORMAT(//6X,'Routine ETPOTS, version 1996.08.05.'/ 1 6x,'Tidal waves from tidal potential catalogue.'/ 2 6X,A20,' tidal potential catalogue is used.'/) 17002 FORMAT(A30,' is existing') 17006 FORMAT(I6,1X,A2,I2,11I3,F12.8,2F12.0,2F10.0,1X,A4) 17008 FORMAT(1X,I4,10I2,4F7.5,F8.4,F9.4,F12.8,1X,A4/F7.5,F8.6,F9.6) 17010 FORMAT(//6x,' Number of waves read from file is :',I6/ 1 6x,' Number of waves above limit is :',I6/ 1 6x,' Number of waves to be used is :',I6/) 17011 FORMAT(I5,1X,A2,2I3,3F10.5,2X,A6,2X,F10.6) 17028 FORMAT(8A10) 17029 FORMAT(6X,8A10) 17030 FORMAT(///6x,'***** Routine ETPOTS finished execution.'/) 17050 FORMAT(/ 1 6x,'***** Error in routine ETPOTS.'/ 2 6x,'***** The current number of waves:',I5,' exceeds the ', 3 'maximum number of waves:',I5/ 4 6x,'***** Routine ETPOTS stops the execution.'/) END C SUBROUTINE GEOEXT(IUN16,IRESET,DEXTIM,DEXTOT) C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C ! C Routine GEOEXT, version 1996.08.05 Fortran 77/90. ! C ! C === MS-DOS version for LAHEY-compiler =================== ! C ! C The routine GEOEXT computes the actual job time and writes ! C the actual execution time on printer output unit IUN6. ! C For the first call of routine GEOEXT, the actual jobtime will ! C be computed (in secs since midnight) and stored. For the next ! C call(s) of routine GEOEXT, the actual jobtime will be computed ! C and the execution time (actual jobtime minus jobtime of the ! C first call of routine GEOEXT) will be printed. ! C ! C Input parameter description: ! C ---------------------------- ! C ! C IUN16: formatted printer unit. ! C IRESET: DEXTIM will be resetted, if IRESET=1. ! C ! C Output parameter description: ! C ----------------------------- ! C ! C DEXTIM: actual jobtime in seconds (time elapsed from the ! C last call of routine GEOEXT with IRESET=1 to the ! C actual call of routine GEOEXT), double precision. ! C DEXTOT: total jobtime in seconds (time elapsed from the ! C first call of routine GEOEXT), double precision. ! C ! C Used routines: ! C -------------- ! C ! C SYSTEM-CLOCK ! C ! C Program creation: 1979.08.30 by Hans-Georg Wenzel, ! C Black Forest Observatory, ! C Universitaet Karlsruhe, ! C Englerstr. 7, ! C D-76128 KARLSRUHE, ! C Germany. ! C Tel.: 0721-6082301. ! C FAX: 0721-694552. ! C e-mail: wenzel@gik.bau-verm.uni-karlsruhe.de ! C Last Modification: 1996.08.05 by Hans-Georg Wenzel. ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IMPLICIT DOUBLE PRECISION (D) C MSFOR: INTEGER*2 IH,IM,IS,IS100 DATA IFIRST/1/ SAVE DTIME1 IF(IRESET.NE.1) GOTO 6003 C MSFOR: CALL GETTIM(IH,IM,IS,IS100) C MSFOR: DTIME1=DBLE(IS+IM*60+IH*3600)+0.01*FLOAT(IS100) C LAHEY 90: CALL SYSTEM_CLOCK(IC,ICR) DTIME1=DBLE(IC)/DBLE(ICR) C UNIX: DTIME1=DBLE(SECNDS(RDUMMY)) WRITE(IUN16,17001) DEXTIM=0.D0 DEXTOT=0.D0 IF(IFIRST.EQ.1) THEN DTIME0=DTIME1 IFIRST=0 ENDIF IRESET=0 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 SYSTEM_CLOCK(IC,ICR) DTIME2=DBLE(IC)/DBLE(ICR) C UNIX: DTIME2=DBLE(SECNDS(RDUMMY)) DEXTIM=DTIME2-DTIME1 DEXTOT=DTIME2-DTIME0 WRITE(IUN16,17002) DEXTIM C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C Format statements: ! C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 17001 FORMAT(6x,'First call of routine GEOEXT, version 1996.08.05.') 17002 FORMAT(/6x,'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