      subroutine admfnd(amp,f,p,n)
c  Returns the amplitude, frequency, and phase of a set of tidal
c  constituents. n is input as the number wanted, and returned as the number
c  provided.  The constituents used are stored in the arrays idd (doodson
c  number) and tamp (Cartwright-Edden amplitude).  The actual amp and
c  phase of each of these are determined by spline interpolation of the
c  real and imaginary part of the admittance, as specified at a subset
c  of the constituents.
c  The phase is determined for a time set in common block /date/ (see
c  subroutine tdfrph), outside of this subroutine.
c   Note that the arrays f and p must be specified as double precision.
c
c   $$$$calls tdfrph,spline,eval,shells
c
      character*1 ipht
      double precision f,p,fr,pr
      dimension amp(1),f(1),p(1)
c  arrays containing information about all stored constituents
      dimension idd(6,67),tamp(67)
c  arrays containing information about the subset whose amp and phase may
c  be specified, and scratch arrays for the spline routines
c   at most 20 constituents may be specified.
      dimension rl(20),aim(20),rf(20),key(20),scr(20),
     1 zdi(20),zdr(20),di(20),dr(20),sdi(20),sdr(20),tdi(20),tdr(20),
     2 idtst(6)
      data dtr/.017453293/,pi/3.14159265/
      data nt/67/,ncon/20/
      data rl/20*0.0/,aim/20*0.0/,rf/20*0.0/
      data zdi/20*0.0/,zdr/20*0.0/,di/20*0.0/,dr/20*0.0/
      data sdi/20*0.0/,sdr/20*0.0/,tdi/20*0.0/,tdr/20*0.0/
      data tamp/
     1  .63192, .29400, .12099, .07966, .02383,-.02358, .02298, .01932,
     2 -.01786, .01720, .01601, .00467,-.00466,-.00451, .00447, .00447,
     3  .00259,-.00246,-.00218, .00197, .00195, .00192,-.00190, .00180,
     4  .36878,-.26221,-.12203,-.05020, .05001,-.04945, .02062, .02062,
     5  .01129,-.00954,-.00947,-.00802, .00741,-.00730, .00723,-.00714,
     6 -.00664, .00525, .00414, .00409, .00395, .00394, .00343, .00342,
     7  .00293, .00289, .00216, .00194,-.00194,-.00180,-.06663,-.03518,
     8 -.02762,-.01276,-.00673,-.00583,-.00529,-.00288,-.00258, .00765,
     9  .00210, .00100,-.00043/
      data idd/
     1 2, 0, 0, 0, 0, 0,  2, 2,-2, 0, 0, 0,  2,-1, 0, 1, 0, 0,  2, 2, 0,
     2 0, 0, 0,  2, 2, 0, 0, 1, 0,  2, 0, 0, 0,-1, 0,  2,-1, 2,-1, 0, 0,
     3 2,-2, 2, 0, 0, 0,  2, 1, 0,-1, 0, 0,  2, 2,-3, 0, 0, 1,  2,-2, 0,
     4 2, 0, 0,  2,-3, 2, 1, 0, 0,  2, 1,-2, 1, 0, 0,  2,-1, 0, 1,-1, 0,
     5 2, 3, 0,-1, 0, 0,  2, 1, 0, 1, 0, 0,  2, 2, 0, 0, 2, 0,  2, 2,-1,
     6 0, 0,-1,  2, 0,-1, 0, 0, 1,  2, 1, 0, 1, 1, 0,  2, 3, 0,-1, 1, 0,
     7 2, 0, 1, 0, 0,-1,  2, 0,-2, 2, 0, 0,  2,-3, 0, 3, 0, 0,  1, 1, 0,
     8 0, 0, 0,  1,-1, 0, 0, 0, 0,  1, 1,-2, 0, 0, 0,  1,-2, 0, 1, 0, 0,
     9 1, 1, 0, 0, 1, 0,  1,-1, 0, 0,-1, 0,  1, 2, 0,-1, 0, 0,  1, 0, 0,
     1 1, 0, 0,  1, 3, 0, 0, 0, 0,  1,-2, 2,-1, 0, 0,  1,-2, 0, 1,-1, 0,
     2 1,-3, 2, 0, 0, 0,  1, 0, 0,-1, 0, 0,  1, 1, 0, 0,-1, 0,  1, 3, 0,
     3 0, 1, 0,  1, 1,-3, 0, 0, 1,  1,-3, 0, 2, 0, 0,  1, 1, 2, 0, 0, 0,
     4 1, 0, 0, 1, 1, 0,  1, 2, 0,-1, 1, 0,  1, 2,-2, 1, 0, 0,  1, 0, 2,
     5-1, 0, 0,  1,-1, 2, 0, 0, 0,  1, 3,-2, 0, 0, 0,  1, 1, 1, 0, 0,-1,
     6 1, 1,-1, 0, 0, 1,  1, 4, 0,-1, 0, 0,  1, 0,-2, 1, 0, 0,  1,-4, 2,
     7 1, 0, 0,  1,-2, 2,-1,-1, 0,  0, 2, 0, 0, 0, 0,  0, 1, 0,-1, 0, 0,
     8 0, 2, 0, 0, 1, 0,  0, 3, 0,-1, 0, 0,  0, 1,-2, 1, 0, 0,  0, 2,-2,
     9 0, 0, 0,  0, 3, 0,-1, 1, 0,  0, 2, 0,-2, 0, 0,  0, 2, 0, 0, 2, 0,
     1 3,0,0,0,0,0,  3,-1,0,1,0,0,  3,2,0,0,0,0,  3,0,0,0,-1,0/
      k = 0
      nlp = 0
      ndi = 0
      nsd = 0
      ntd = 0
      read(5,104) ipht
 104  format(a1)
      if(ipht.eq.'l') read(5,106) rlong
 106  format(f8.0)
 13   read(5,108,end=17) idtst,ramp,rph
 108  format(6i2,2f9.0)
c  stop reading on a flag or if arrays are full
      if(idtst(1).eq.-1.or.k.eq.ncon) go to 17
c  see if Doodson numbers match
      do 15 kk=1,nt
      ii = 0
      do 14 i=1,6
 14   ii = ii + iabs(idd(i,kk)-idtst(i))
      if(ii.eq.0) go to 16
 15   continue
c  have a match - put line into array
 16   if(ii.eq.0.and.k.lt.ncon) then
         k = k + 1
         rl(k) = ramp*cos(dtr*rph)/abs(tamp(kk))
         aim(k)= ramp*sin(dtr*rph)/abs(tamp(kk))
c  Now have real and imaginary parts of admittance, scaled by C-E
c  amplitude. Admittance phase is whatever was used in the original
c  expression. (Usually phase is given relative to some reference
c  but amplitude is in absolute units). Next get frequency.
         call tdfrph(idd(1,kk),fr,pr)
         rf(k) = fr
      endif
      go to 13
c  done reading in constituents
c  have specified admittance at a number of points. sort these by frequency
c  and separate diurnal and semidiurnal, recopying admittances to get them in
c  order
 17   call shells(rf,key,k)
      do 21 i=1,k
      if(rf(i).lt.0.5) nlp = nlp + 1
      if(rf(i).lt.1.5.and.rf(i).gt.0.5) ndi = ndi + 1
      if(rf(i).lt.2.5.and.rf(i).gt.1.5) nsd = nsd + 1
      if(rf(i).lt.3.5.and.rf(i).gt.2.5) ntd = ntd + 1
 21   scr(i) = rl(key(i))
      do 23 i=1,k
      rl(i) = scr(i)
 23   scr(i) = aim(key(i))
      do 25 i=1,k
 25   aim(i) = scr(i)
c  now set up splines (8 cases - four species, each real and imag)
c  we have to allow for the case when there are no constituent amplitudes
c  for the long-period or terdiurnal.
      if(nlp.ne.0) call spline(nlp,rf,rl,zdr,scr)
      if(nlp.ne.0) call spline(nlp,rf,aim,zdi,scr)
      call spline(ndi,rf(nlp+1),rl(nlp+1),dr,scr)
      call spline(ndi,rf(nlp+1),aim(nlp+1),di,scr)
      call spline(nsd,rf(nlp+ndi+1),rl(nlp+ndi+1),sdr,scr)
      call spline(nsd,rf(nlp+ndi+1),aim(nlp+ndi+1),sdi,scr)
      if(ntd.ne.0) 
     1   call spline(ntd,rf(nlp+ndi+nsd+1),rl(nlp+ndi+nsd+1),tdr,scr)
      if(ntd.ne.0) 
     1   call spline(ntd,rf(nlp+ndi+nsd+1),aim(nlp+ndi+nsd+1),tdi,scr)
c  evaluate all harmonics using the interpolated admittance
      n = min0(n,nt)
      j = 1
      do 31 i=1,n
      if((idd(1,i).eq.0.and.nlp.eq.0).or.(idd(1,i).eq.3.and.ntd.eq.0))
     1  go to 31
      call tdfrph(idd(1,i),f(j),p(j))
c  equilibrium phase corrections (including to local phase)
      if(idd(1,i).eq.1) p(j) = p(j) + 90.
      if(idd(1,i).eq.3) p(j) = p(j) + 180.
      if(ipht.eq.'l') p(j) = p(j) + idd(1,i)*rlong
      sf = f(j)
      if(idd(1,i).eq.0) re = eval(sf,nlp,rf,rl,zdr)
      if(idd(1,i).eq.0) am = eval(sf,nlp,rf,aim,zdi)
      if(idd(1,i).eq.1) re = eval(sf,ndi,rf(nlp+1),rl(nlp+1),dr)
      if(idd(1,i).eq.1) am = eval(sf,ndi,rf(nlp+1),aim(nlp+1),di)
      if(idd(1,i).eq.2) re =
     1   eval(sf,nsd,rf(nlp+ndi+1),rl(nlp+ndi+1),sdr)
      if(idd(1,i).eq.2) am =
     1   eval(sf,nsd,rf(nlp+ndi+1),aim(nlp+ndi+1),sdi)
      if(idd(1,i).eq.3) re =
     1   eval(sf,ntd,rf(nlp+ndi+nsd+1),rl(nlp+ndi+nsd+1),tdr)
      if(idd(1,i).eq.3) am =
     1   eval(sf,ntd,rf(nlp+ndi+nsd+1),aim(nlp+ndi+nsd+1),tdi)
      amp(j) = tamp(i)*sqrt(re**2+am**2)
      p(j) = p(j) + atan2(am,re)/dtr
      j = j + 1
 31   continue
      n = j - 1
      return
      end
