      subroutine tdfrph(idood,freq,phase)
      double precision freq,phase,fr,ph,tm,rft
      dimension fr(6,2),ph(6,2),idood(6),rft(2)
c
c   Given the Doodson number of a tidal constituent (in idood), returns
c  the frequency and phase.  Phase is returned in degrees and frequency
c  in cycles/day.
c
c   Note that phases must be decreased by 90 degrees if the sum of the order
c  and the species number is odd (as for the 2nd degree diurnals, and 3rd
c  degree low-frequency and semidiurnals).
c   These phases may need further adjustment to allow for the spherical
c  harmonic normalization used; e.g. for that used for the potential by
c  Cartwright and Tayler, 180 degrees must be added for (species,order)
c  = (1,2), (1,3), or (3,3).
c
c  Common block date stores time information. phase at given time is computed
c  from frequency, and phase at reference times (1924, day 142.5, or 1960,
c  day 143.5). Times before 1900 will cause a stop.
c  Frequencies and phases are from Cartwright and Tayler, Geophys. J.,
c  vol. 23, p 56 (1971)
c
c       calls only system routines
c
      common/date/iyr,iday,ihr,imin,isec
      data rft/8906.5,22056.5/,
     2 fr/1.d0,3.66011025d-2,2.7379092d-3,3.094554d-4,1.470941d-4,
     3    1.3708d-7,
     4    1.d0,3.66011023d-2,2.7379092d-3,3.094548d-4,1.47094d-4,
     5    1.3708d-7/,
     6 ph/1.8d2,272.60245d0,58.85684d0,246.60455d0,212.47704d0,
     7    281.64011d0,
     8    1.8d2, 22.22101d0,60.11923d0,271.56503d0,188.82048d0,
     9    282.25919d0/
c  find the time in days from 1 jan 1900 (if before 1900 stop)
      if(iyr.lt.1900) stop
      tm = iday - 1 + 365.d0*(iyr-1900) + (iyr-1901)/4
c  see if we want to use the first or second set of numbers
      ind = tm/15482. + 1.
      if(ind.gt.2) ind = 2
      fday = ihr/24.d0 + imin/1440.d0 + isec/86400.d0
      tm = tm + fday - rft(ind)
      freq = idood(1)*(fr(1,ind) - fr(2,ind) + fr(3,ind))
      phase = idood(1)*(ph(1,ind) - ph(2,ind) + ph(3,ind))
      do 7 i = 2,6
      freq = freq + fr(i,ind)*idood(i)
 7    phase = phase + ph(i,ind)*idood(i)
      phase = phase + 360.d0*freq*tm
c adjust phases to fall in the positive range 0 to 360
      phase = dmod(phase,360.d0)
      if(phase.lt.0.d0) phase = phase + 360.d0
      return
      end
