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