c  program hartid
c    finds tides by summing harmonics; the amp and phase of these are
c  specified in subroutine admfnd
c
c    $calls admget,tdfrph,recurs,spline,shells
c
      character*40 dumm
      double precision f,p,dr,pi,scr
      dimension a(67),f(67),p(67),x(600),hc(134),scr(201),wf(67)
      common/date/it(5)
      data nt/67/,dr/.01745329252d0/,pi/3.1415926535898d0/,irli/1/
      data luo/6/
      if(iargc().lt.1.or.iargc().gt.7) then
	 write(luo,100)
 100  format('Usage: hartid y d h m s num samp',/,
     1  '  harmonics file read from standard input',/,
     2  '  results written to standard output')
	 stop
      endif
      call getarg(1,dumm)
      read(dumm,102) it(1)
 102  format(i4)
      call getarg(2,dumm)
      read(dumm,102) it(2)
      call getarg(3,dumm)
      read(dumm,102) it(3)
      call getarg(4,dumm)
      read(dumm,102) it(4)
      call getarg(5,dumm)
      read(dumm,102) it(5)
      call getarg(6,dumm)
      read(dumm,104) irnt
 104  format(i6)
      call getarg(7,dumm)
      read(dumm,106) samp
 106  format(f7.0)
      call admfnd(a,f,p,nt)
c  set up for first recursion, and normalize frequencies
      do 23 i=1,nt
      p(i) = dr*p(i)
      f(i) = samp*pi*f(i)/43200.d0
 23   wf(i) = f(i)
 31   irhi = min(irli+599,irnt)
      np = irhi - irli + 1
c set up harmonic coefficients, compute tide, and write out
      do 33 i=1,nt
      hc(2*i-1) = a(i)*dcos(p(i))
 33   hc(2*i)  = -a(i)*dsin(p(i))
      call recurs(x,np,hc,nt,wf,scr)
      write(6,110) (x(i),i=1,np)
 110  format(6g13.5)
      if(irhi.eq.irnt) go to 37
      irli = irhi + 1
c  reset phases to the start of the new section
      do 36 i=1,nt
      p(i) = p(i) + np*f(i)
 36   p(i) = dmod(p(i),2.d0*pi)
      go to 31
 37   continue
      stop
      end
