        program tidiff
c       tide signal zur Seismographeneichung, U. Noell April 1993
c   Gibt die zeitliche Ableitung des tide signals aus, entsprechend
c  der Antwort eines VBB - Seismographen
c  this program uses TIDEPT to compute a single component of the earth
c  tide signal on an elastic earth.
c  besides vital variables,format statements 108,109 should be modified
      dimension xn(120000),xe(120000),xnv(120000),xev(120000)
      dimension xz(120000),xzv(120000)
      double precision ts,dt,t,tim69
      integer unit,year,month,day
      character prolog(5)*50,text*50
      logical ex
      inquire(file='tides.par',exist=ex)
      if(ex.eq..true.) then
      open(1,file='tides.par')
      unit=1
      else
      unit=5
      endif
      open(unit=2,file='tidiff.out')
      write(*,*)
      write(*,*) '>> If file tides.par exists, then input data are read'
      write(*,*) '>> from it.  Otherwise, input is interactive.'
      write(*,*) '>> Input can also be piped in with tides < inputfile'
      write(*,*)
      write(*,*) 'latitude?'
      read(unit,*) glat
      write(*,*) 'latitude=',glat
      write(2,*) 'latitude=',glat
      write(*,*) 'longitude?'
      read(unit,*) glong
      write(*,*) 'longitude=',glong
      write(2,*) 'longitude=',glong
      write(*,*) 'year,month,day (GMT)? '
      read(unit,*) year,month,day
      write(*,*) 'year,month,day (GMT) =',year,month,day
      write(2,*) 'year,month,day (GMT) =',year,month,day      
      write(*,*) 'how many days? (sampled every quarter hour):'
      read(unit,*) ntage
      dt=0.25d0
      ndata=min(120000,nint(24./dt*ntage))
      write(*,*) ntage," days, ",ndata," data points. dt=15 minutes."
      write(2,*) ntage," days, ",ndata," data points. dt=15 minutes."
      write(prolog(4),*) 
     &    ntage," days, ",ndata," data points. dt=15 minutes."
      write(*,*) 
     &  "time in hours GMT, tidal derivative in microgals per hour"
        write(2,*) 
     &  "time in hours GMT, tidal derivative in microgals per hour"
      write(prolog(5),*) 
     &  "time in hours GMT, tidal derivative in microgals per hour"
      write(2,*) 
     &'     time ,  z-component , n-component  ,  e-component ' 
C      glat= 48.7708
C      glong= 9.19333
c  east longitude in degrees
      azim=90.000
c  azimuth in degrees
C      wink=2.74017
c      strainmeter:east from north
c      tiltmeter:north from east,i.e. 0.0 degrees is east,90 degrees is
c              north tilt
C      opt=2
C      opt=3
C            opt=1
      lorder=0
c  gravity=1,tilt=2,strain=3
c      dt=1.00d0
c  timestep in hours
      ts=tim69(year,month,day)
C      ts= tim69(1991,7,25)
c  first point,time in hours from jan 01,1969,0.0hrs GMT
c     scale= 1.00e+06
c      scale= 0.206264806e+09
      scale= 1.00e+06
c   if scale =1.0e+6 then gravity in microgals,tilt in microrads,strain
c       in microstrain.for nanorads and nanostrains resp. use 1.0e+09
c       for tilt in msecs arc use  0.206264806e+09
C      ndata = 32
      t=ts
      do  1  i=1,ndata
      call tidept(t,xz(i),glat,glong,azim,scale,lorder,1)
      call tidept(t,xn(i),glat,glong,azim,scale,lorder,2)
      azm2 = 0.00
      call tidept(t,xe(i),glat,glong,azm2,scale,lorder,2)
C      azm2 = 45.00
C      call tidept(t,str1,glat,glong,azm2,scale,lorder,opt)
C      azm2 = -45.00
C      call tidept(t,str2,glat,glong,azm2,scale,lorder,opt)
C      shear=(str1-str2)/2.00
C      dil=(x(i)+str)*(0.66666)
C      senk=str*cos(wink)**2+x(i)*sin(wink)**2-shear*sin(2.*wink)
C      write(8,112) i,dil,senk
      t=t+dt
   1  continue
      do 11 i=2,ndata-1
      xzv(i)=-1.*(xz(i+1)-xz(i-1))/2./dt
      xnv(i)=-981.*(xn(i+1)-xn(i-1))/2./dt
      xev(i)=-981.*(xe(i+1)-xe(i-1))/2./dt
      write(2,112) i*dt,xzv(i),xnv(i),xev(i)
11    continue
      close(2)
      si=900.
      tmin=15.
      tsec=0.
      text='tidal derivative EAST from program TIDIFF         '
      call outform('tidiff.e',text,prolog,ndata,si,tmin,tsec,xev)
      text='tidal derivative NORTH from program TIDIFF        '
      call outform('tidiff.n',text,prolog,ndata,si,tmin,tsec,xnv)
      text='tidal derivative VERTICAL from program TIDIFF     '
      call outform('tidiff.v',text,prolog,ndata,si,tmin,tsec,xzv)

900      format('  i    east-strain  north-strain  shear-strain ')
901      format('  i    nord-tilt    east-tilt ')
902      format('  i    Dilatation    strain senkrecht auf Kluft ')
 112  format(4(f12.4,2x))
      
c  generate plot-parameter file
      open(4,file='winplot.par')
      write(4,*) '20, 3, 24, 15, 1,',ndata, ', 0.8'
      write(4,'(a)') 'tidiff.e'
      write(4,'(a)') 'tidiff.n'
      write(4,'(a)') 'tidiff.v'
      write(4,'(a)') 'synthetic tidal derivative from program TIDIFF'
      write(*,*) 'plot-parameter file winplot.par was generated'
 
      write(*,*) 
     & 'The three time series were also written into file tidiff.out'
      stop
      end

      subroutine outform(name,text,prolog,n,dt,tmin,tsec,x)
c  formatted output
      dimension x(120000)
      character name*8,text*50,iform*20,prolog(5)*50
      xmax=0.
      do 2 j=1,n
    2 xmax=max(xmax,abs(x(j)))
      nvor=int(log10(max(xmax,2.))+1.)
      ndec=max(0,10-nvor)
      write(iform,3) ndec
    3 format('(5f13.',i1,')')
      write(6,5) n,trim(name)
    5 format(' writing ',i7,' samples into file ',a8,' in SEIFE format')
      open(8,file=trim(name))
      write(8,'(a)') text
      do i=1,5
        write(8,'(a5,a50)') '%    ',prolog(i)
      enddo
      write(8,1) n,iform,dt,tmin,tsec
    1 format(i10,a20,f10.3,2f10.3)
      write(8,iform) (x(j),j=1,n)
      close(8)
      return
      end


        Double Precision Function tim69(myr,mon,mday)
c  function to return the number of days from the beginning of 1969
c  to a given later date.
c  intended to be used to find the input variable for the program TIDEPT
c  and its relatives, which need time in hours from January 1,1969,
c  0.0 hours G.M.T.
c  this version goes wrong after 2100 and before 1800
c
      double precision st69,hours
      dimension iam(12)
      data iam/0,31,59,90,120,151,181,212,243,273,304,334/
      lday=mday-1
      kyr=myr-1969
      lyr=kyr
      if(lyr.lt.0)  lyr=lyr-3
      st69=kyr*365+iam(mon)+lday+lyr/4
      if(myr.eq.myr/4*4.and.mon.gt.2) st69=st69+1.0d0
      if(myr.le.1900.and.mon.le.2)  st69=st69+1.0d0
      hours = st69*24.00d0
      tim69=hours
      return
      end


      SUBROUTINE TIDEPT (TIME,RESULT,GLAT,GLONG,ORIENT,SCALE,
     * LORDER,OPT)
C    BERGER,EVANS,MCKENZIE
C MODIFIED BY YOUNG 16.1.78 TO RETURN STRAIN VALUE AT ONE TIME
C ONLY.  SEE SUBROUTINE EARTHTD FOR DETAILS OF VARIABLES.
C     STRAINMETER AZIMUTH IS POSITIVE EAST FROM NORTH
      REAL LAMDA,   LL,         H,K,L,KR,MU,MUT,MUL,MUTT,MULL,MUTL,ELL,
     1  ETT,ETL
      DOUBLE PRECISION  SAW,CAW,RBARR
      DOUBLE PRECISION CZ,SZ,Z,X,THFPSI,CNM,SNM,CW,SW,CNU,SNU,NU
      DOUBLE PRECISION TIME,TT,PSIG,PSIM,HR,AW
      DOUBLE PRECISION  PI20,HALFPI,TWOPI,AZT,HS,T2,PS,T,PSIS,HM,PM,NM,
     1  ES,SIG,LS,LM,TANP
      DOUBLE PRECISION  DATANZ
      INTEGER OPT
      DIMENSION  H(6),K(6),L(6),KR(7),P(7),PP(7),PPP(7)
      H(1)=0.612
      H(2)=0.290
      H(3)=0.175
      H(4)=0.129
      H(5)=0.107
      H(6)=0.095
      K(1)=0.302
      K(2)=0.093
      K(3)=0.042
      K(4)=0.025
      K(5)=0.017
      K(6)=0.013
      L(1)=0.083
      L(2)=0.014
      L(3)=0.010
      L(4)=0.008
      L(5)=0.007
      L(6)=0.005
      JS=2
      JE=3
      IF(LORDER.EQ.0)  GO TO5
      JE=LORDER
      JS=JE
  5   CONTINUE
      GTSUN=0.
      TTSUN=0.
      STNSUN=0.
      LAMDA=GLONG
      THETA=90.0-GLAT
      AZS=-ORIENT
      AZT=-ORIENT
      PI20=62.83185307169D0
      HALFPI=1.570796326793D0
      TWOPI=6.283185307169D0
      RE=6.37122E+8
      G=982.02
      AZT=AZT*0.0174532925199D0
      AZS=AZS*0.0174532925199D0
      CAZT=COS(AZT)
      SAZT=SIN(AZT)
      CAZS=COS(AZS)
      SAZS=SIN(AZS)
      TT=TIME+604848.0D0
      THETA=THETA*0.0174532925199D0
      LAMDA=LAMDA*0.0174532925199D0
      CTHETA=COS(THETA)
      STHETA=SIN(THETA)
      HR=DMOD(TIME,24.0D0)
      T=(TT+12.0D0)/876600.0D0
      T2=T*T
      HS=4.881627982482D0+628.3319508731D0*T
      HS=HS+0.523598775578D-5*T2
      HS=DMOD(DMOD(HS,PI20)+PI20,PI20)
      PS=4.908229466993D0+T*0.0300052641669D0
      PS=PS+0.790246300201D-5*T2
      ES=0.01675104D0- 0.00004180D0*T-0.000000126D0*T2
      PSIG=0.2617993877971D0*(HR-12.0D0)+HS
      IF(LORDER.NE.0.AND.LORDER.NE.2)   GO TO 82
      LS=HS+2.0D0*ES*DSIN(HS-PS)+1.25D0*ES**2*DSIN(2.0D0*(HS-PS))
      CZ=0.397980654630D0*DSIN(LS)
      SZ=DSQRT(1.0D0-CZ**2)
      Z=DATANZ(SZ,CZ)
      X=0.5D0*(LS+Z-HALFPI)
      THFPSI=1.52386101015D0*(DSIN(X)/DCOS(X))
      PSIS=2.0D0*DATAN(THFPSI)
      PSIS=DMOD(PSIS+TWOPI,TWOPI)
      RBARR=1.D0+ES*DCOS(HS-PS)+ES**2*DCOS(2.D0*(HS-PS))
      LL=PSIS-PSIG
      SLL=SIN(LAMDA-LL)
      CLL=COS(LAMDA-LL)
      CALFA=CTHETA*CZ+STHETA*SZ*CLL
      XI=4.2635E-5*RBARR
      CC=2.120E14*XI
      KR(2)=CC*XI*XI
      KR(3)=KR(2)*XI
      KR(4)=KR(3)*XI
      MU=CALFA
      P(2)=0.5*(3.*MU*MU-1.)
      P(3)=0.5*MU*(5.*MU*MU-3.)
      P(4)=0.25*(7.*MU*P(3)-3.*P(2))
      PP(2)=3.*MU
      PP(3)=1.5*(5.*MU*MU-1.)
      PP(4)=0.25*(7.*(P(3)+MU*PP(3))-3.*PP(2))
      PPP(2)=3.
      PPP(3)=15.*MU
      PPP(4)=7.5*(7.*MU*MU-1.)
      MUT=-STHETA*CZ+CTHETA*SZ*CLL
      MUTT=-MU
      MUL=-STHETA*SZ*SLL
      MULL=-STHETA*SZ*CLL
      MUTL=-CTHETA*SZ*SLL
      J=2
      GO TO(50,60,70),OPT
 50   GTSUN=0.
      DEL=1.+2.*H(J-1)/J-(J+1)*K(J-1)/J
      GTSUN=GTSUN+DEL*J*KR(J)*P(J)*G/RE
      GO TO 80
 60   TTSUN=0.
      DIM=1.+K(J-1)-H(J-1)
      TTSUN=TTSUN+(DIM*KR(J)*PP(J)/RE)*(MUT*SAZT+MUL*CAZT/STHETA)
      GO TO 80
 70    ELL=0.
      ETL=0.
      ETT=0.
      M=J-1
      ELL=ELL+KR(J)*(L(M)*(PP(J)*MUT*CTHETA/STHETA+(PPP(J)*
     1MUL*MUL+PP(J)*MULL)/(STHETA*STHETA))+H(M)*P(J))
      ETT=ETT+KR(J)*(L(M)*(PPP(J)*MUT*MUT+PP(J)*MUTT)+
     1  H(M)*P(J))
      ETL=ETL+2.*(KR(J)*L(M)*(PPP(J)*MUL*MUT+PP(J)*MUTL-
     1  CTHETA*PP(J)*MUL/STHETA))/STHETA
      STNSUN=(ELL*SAZS*SAZS+ETL*SAZS*CAZS+ETT*CAZS*CAZS)/RE
 80   CONTINUE
 82   HM=4.720008893870D0+8399.709274530D0*T
      HM=HM+0.345575191895D-4*T2
      HM=DMOD(DMOD(HM,PI20)+PI20,PI20)
      PM=5.835152597865D0+71.01804120839D0*T
      PM=PM-0.180108282526D-3*T2
      NM=4.523601611611D0-33.75714629423D0*T
      NM=NM+0.362640633469D-4*T2
      CNM=DCOS(NM)
      SNM=DSIN(NM)
      CW=0.91369D0-0.03569D0*CNM
      SW=DSQRT(1.0D0-CW**2)
      SNU=0.08968D0*SNM/SW
      CNU=DSQRT(1.0D0-SNU**2)
      NU=DATAN(SNU/CNU)
      SAW=0.39798D0*SNM/SW
      CAW=CNM*CNU+0.91739D0*SNM*SNU
      AW=2.0D0*DATAN(SAW/(1.0D0+CAW))
      AW=DMOD(AW+TWOPI,TWOPI)
      SIG=HM-NM+AW
      LM=SIG+0.109801D0*DSIN(HM-PM)+0.003768D0*DSIN(2.0D0*(HM-PM))
     1 -0.224412D0*ES*DSIN(HS-PS)+0.02045D0*DSIN(HM-2.0D0*HS+PM)
     2 +0.010809D0*DSIN(2.0D0*(HM-HS))+0.000653D0*DSIN(3.0D0*HM-2.0D0*HS
     3  -PM) +0.000451D0*DSIN(2.0D0*HM-3.0D0*HS+PS)
      CZ=DSIN(LM)*SW
      WM=DATANZ(SW,CW)
      SZ=DSQRT(1.0D0-CZ**2)
      ZM=DATANZ(SZ,CZ)
      X=0.5D0*(LM+ZM-HALFPI)
      TANP=(DSIN(X)/DCOS(X))*DSIN(0.5D0*(HALFPI+WM))/DSIN(0.5D0*(HALFPI-
     1  WM))
      PSIM=2.0D0*DATAN(TANP)
      PSIM=DMOD(PSIM+TWOPI,TWOPI)
      RBARR=1.0D0+0.05467D0*DCOS(HM-PM)+0.003011D0*DCOS(2.0D0*(HM-PM))
     1 -0.000141D0*DCOS(HS-PS)+0.009271D0*DCOS(HM-2.0D0*HS+PM)
     2+0.007759D0*DCOS(2.0D0*(HM-HS))+0.000633D0*DCOS(3.D0*HM-2.D0*HS-PM
     3) +0.000328D0*DCOS(2.D0*HM-3.D0*HS+PS)
      LL=NU+PSIM-PSIG
      SLL=SIN(LAMDA-LL)
      CLL=COS(LAMDA-LL)
      CALFA=CTHETA*CZ+STHETA*SZ*CLL
      XI=1.65933E-2*RBARR
      CC=7.834E6*XI
      KR(2)=CC*XI*XI
      J=3
      KR(J)=KR(J-1)*XI
      MU=CALFA
      P(1)=MU
      P(2)=0.5*(3.*MU*MU-1.)
      J=2
      P(J+1)=((2*J+1)*MU*P(J)-J*P(J-1))/(J+1)
      PP(1)=1.
      PP(2)=3*MU
      PPP(1)=0.
      PPP(2)=3.
      J=2
      PP(J+1)=((2*J+1)*(P(J)+MU*PP(J))-J*PP(J-1))/(J+1)
      PPP(J+1)=((2*J+1)*(2*PP(J)+MU*PPP(J))-J*PPP(J-1))/(J+1)
      MUT=-STHETA*CZ+CTHETA*SZ*CLL
      MUTT=-MU
      MUL=-STHETA*SZ*SLL
      MULL=-STHETA*SZ*CLL
      MUTL=-CTHETA*SZ*SLL
      GO TO (100,110,120),OPT
 100   GTMUN=0.
      DO  105  J=JS,JE
      DEL=1.+2*H(J-1)/J-(J+1)*K(J-1)/J
 105  GTMUN=GTMUN+DEL*J*KR(J)*P(J)*G/RE
      GO TO 130
 110   TTMUN=0.
      DO  115  J=JS,JE
      DIM=1.+K(J-1)-H(J-1)
 115  TTMUN=TTMUN+(DIM*KR(J)*PP(J)/RE)*(MUT*SAZT+MUL*CAZT/STHETA)
      GO TO 130
 120  ELL=0.
      ETL=0.
      ETT=0.
      DO 125 J=JS,JE
      M=J-1
      ELL=ELL+KR(J)*(L(M)*(PP(J)*MUT*CTHETA/STHETA+(PPP(J)*
     1  MUL*MUL+PP(J)*MULL)/(STHETA*STHETA))+H(M)*P(J))
      ETT=ETT+KR(J)*(L(M)*(PPP(J)*MUT*MUT+PP(J)*MUTT)+H(M)*P(J))
      ETL=ETL+2*(KR(J)*L(M)*(PPP(J)*MUL*MUT+PP(J)*MUTL-
     1 CTHETA*PP(J)*MUL/STHETA))/STHETA
 125  STNMUN=(ELL*SAZS*SAZS+ETL*SAZS*CAZS+ETT*CAZS*CAZS)/RE
 130  CONTINUE
      GO TO (140,150,160),OPT
140   RESULT=(GTSUN+GTMUN)*SCALE
      GO TO 170
150   RESULT=(TTSUN+TTMUN)*SCALE
      GO TO 170
160   RESULT=(STNSUN+STNMUN)*SCALE
170   CONTINUE
 180  CONTINUE
      RETURN
      END
c  datanz returns values of phase in rads between plus and minus pi
c  has been tested in all quadrants
c  if both y and x are equal to zero,datanz returns pi/2
      double precision function datanz(y,x)
      double precision x,y,halfpi
      halfpi= 2.00d0*datan(1.000d0)
      if(x.ne.0.0d0)  go to 1
      if(y.ge.0.00d0)  datanz=halfpi
      if(y.lt.0.00d0)  datanz=-halfpi
      return
   1  datanz=datan(y/x)
      if(x.lt.0.0d0)  datanz=datanz+2.00d0*halfpi
      if(y.lt.0.0d0.and.x.lt.0.0d0) datanz=datanz-4.00d0*halfpi
      return
      end
