      subroutine tptide( dlat,dlon,time,tide,isdata,u,v,pseudo,
     .                   radiation )
*
* this code is based upon the gstide program used to compute the
*      cartwright&ray geosat tide solution.
*
*
*  name - tptide              name derivation - topex/poseidon derived tide
*
*  function -  to compute the ocean tidal height at a given time
*              and location.
*
*  language - fortran 77
*
*  arguments -
*     name      type  i/o               description
*     ----      ----  ---               -----------
*     dlat       d     i    north latitude (in degrees, -90 to  90)
*
*     dlon       d     i    east longitude (in degrees,   0 to 360).
*
*     time       d     i    desired time, in seconds of mjd.  e.g.,
*                            jan 1, 1986 00:00 is 4011638400.
*
*     tide       d     o    computed tidal height, in cm.
*
*     isdata     l     o    logical denoting whether tide data exist at
*                           desired location.  if false, then tide is
*                           not modified.
*
*     u,v        d     o    the orthoweights interpolated to dlat,dlon
*                           dimensioned (3,2)
*
*     pseudo     l     o    true for c&r style "orthoweights"
*                           this can be passed to admit.f for proper
*                           interpretation of the orthoweights.
*                           set per integer on orthoweight file
*
*     radiation  l     o    true if the correction for radiation
*                           potential is activated.
*                           set per integer on orthoweight file
*
*
*  usage notes -
*     using the user supplied orthoweight file, this routine computes
*     the height of the (ocean + load) tide at the desired location.
*     this quantity is appropriate for use in satellite altimetry.
*     the computed tide is composed of the 30 largest spectral lines
*     within both the diurnal and semidiurnal bands, sufficient for
*     representing all major constituents including nodal modulations.
*
*  processing logic -
*     the tidal height is computed at four grid points surrounding
*     the desired location, with the final result being a bilinear
*     interpolation.
*
*  file references -
*     the orthoweight data are read on initial call to routine:
*     file iunit - global array of orthoweights defining tidal
*                  constants at all valid locations. (see parameter statement)
*
*  important local variables -
*     array uv is loaded with orthoweights as follows:
*         uv(i,j,l,m,n), where
*           i = 1 for u; 2 for v
*           j = 1,2 or 3  (3 complex coeffs per species.)
*           l = tidal species (1 or 2)
*           m = 1 to nx, longitude index number
*           n = 1 to ny, latitude index
*
*     for machines like the cray, the compiler option to ignore double
*         precision should be selected.
*
*  error processing - none.
*
*  technical references -
*     d. cartwright and r. ray, oceanic tides from geosat altimetry,
*        journal of geophysical research, 95, 3069-3090, 1990.
*        (particularly appendix a).
*
*  history -
*   version   date       programmer        change description
*   -------   ----       ----------        ------------------
*     1.0    8/17/90   ray/cartwright  initial version.
*     1.1   12/18/90   ray/cartwright  enhanced documentation;
*                                      use ascii input for portability.
*     1.2   12/31/91   r. ray        changed min w to 0.5 (old value of
*                                    0.2 allowed too much extrapolation)
*     2.0    2/25/94   r. eanes      changed to use t/p derived orthoweights
*     2.1    3/02/94   r. eanes      generalized grid size up to 1 x 1 degree
*                                    and added data statements for the harmonic
*                                    amplitudes and phases (a la the jpl mods)
*                                    interpolate orthoweights instead of tides
*     2.2    3/03/94   r. eanes      return orthoweights (u,v) for other use
*                                    added itype to switch to c&r style orthoweights
*                                    change min w (wmin) back to 0.2.
*     2.3    4/18/94   r. eanes      take out t/p kludge for latitudes > 66.2
*     2.4    5/20/94   r. eanes      change treatment of extrapolating past the southern
*                                    and northernmost grid points.  set wmin to 0.1.
*     2.5                            added r.ray's radiation tide correction
*                                    wmin changed to 0.5 (again)
*     3.0    1/24/95   r. eanes      parameter statements for dimensions and unit
*                                    number.  compact format option for orthoweight file.
*                                    increased dimensions to allow 0.5 x 0.5 deg grid.
*                                    changed wmin to 0.2 to make rsc94b work at gauges.
*

      implicit double precision (a-h,o-z)
      parameter (nxmax=720,nymax=361)
      parameter (iunit=54)
      character*80     title1,title2
      real             uv(2,3,2,nxmax,nymax), undef
      double precision latmin,latmax,lonmin,lonmax
      dimension        u(3,2),v(3,2),iuv(12)
      save
      logical    init,isdata,pseudo,radiation
      data       init/.true./, undef/99999./, eps /1.d-10/, wmin /0.2/


*     on first call, read t/p   -derived ortho-weights
*     ------------------------------------------------
      if (init) then
         init = .false.
         do 4 j=1,nymax
         do 4 i=1,nxmax
            uv(1,1,1,i,j) = undef
    4    continue
         read (iunit,14) title1,title2
         write (6,14) title1,title2
         read (iunit,11) nx,ny,latmin,latmax,lonmin,lonmax,idiff,
     .                   itype,irad,ifmt
c
c...orthoweight file format options
c   idiff - if .eq. 1 then zero tide returned for points outside of domain
c                     and isdata set to true
c   itype - if .eq. 1 then cartwright&ray 1991 pseudo-orthoweights are assumed
c   irad  - if .eq. 1 then the cartwright&ray radiation potential modification is used
c   ifmt  - if .eq. 1 then compact integer format used for orthoweights
c
         if (itype.eq.1) then
            pseudo = .true.  ! c&r style orthoweights
         else
            pseudo = .false. ! corrected orthoweights
         endif
         if (irad.eq.1) then
            radiation = .true.  ! corrected for radiation potential
         else
            radiation = .false. ! not corrected for radiation potential
         endif
         if (nx.gt.nxmax .or. ny.gt.nymax) stop 1   ! must change dimensions
         dx = (lonmax-lonmin)/(nx-1)
         dy = (latmax-latmin)/(ny-1)
         xlats = latmin - dy                    ! southern limit of extrapolation
         xlatn = latmax + dy                    ! northern limit of extrapolation
    8    continue
         if (ifmt.eq.0) then                    ! original c&r format
            read(iunit,12,end=10) i,j,u1,v1,u2,v2,u3,v3,
     .                                u4,v4,u5,v5,u6,v6
            uv(1,1,1,i,j) = u1
            uv(2,1,1,i,j) = v1
            uv(1,2,1,i,j) = u2
            uv(2,2,1,i,j) = v2
            uv(1,3,1,i,j) = u3
            uv(2,3,1,i,j) = v3
            uv(1,1,2,i,j) = u4
            uv(2,1,2,i,j) = v4
            uv(1,2,2,i,j) = u5
            uv(2,2,2,i,j) = v5
            uv(1,3,2,i,j) = u6
            uv(2,3,2,i,j) = v6
         endif
         if (ifmt.eq.1) then                    ! compact integer format
            read(iunit,*,end=10) i,j,iuv
            uv(1,1,1,i,j) =  iuv(1)/1000.d0
            uv(2,1,1,i,j) =  iuv(2)/1000.d0
            uv(1,2,1,i,j) =  iuv(3)/1000.d0
            uv(2,2,1,i,j) =  iuv(4)/1000.d0
            uv(1,3,1,i,j) =  iuv(5)/1000.d0
            uv(2,3,1,i,j) =  iuv(6)/1000.d0
            uv(1,1,2,i,j) =  iuv(7)/1000.d0
            uv(2,1,2,i,j) =  iuv(8)/1000.d0
            uv(1,2,2,i,j) =  iuv(9)/1000.d0
            uv(2,2,2,i,j) = iuv(10)/1000.d0
            uv(1,3,2,i,j) = iuv(11)/1000.d0
            uv(2,3,2,i,j) = iuv(12)/1000.d0
         endif
         go to 8
   10    continue
   11    format(2i4,4f8.3,4i2)
   12    format(2i4,16x,6f8.3,/,24x,6f8.3)
   14    format(a80,/,a80)
      endif
      isdata = .true.
      tide = 0.0d0
      do i=1,3
         do j=1,2
            u(i,j) = 0.d0
            v(i,j) = 0.d0
         enddo
      enddo

*     compute indices for desired position
*     ------------------------------------
      xlat = dlat
      jlat1 = int((xlat - latmin)/dy) + 1
      if (jlat1.eq.0  .and. xlat.gt.xlats) jlat1 = 1        ! allow extrapolation to xlats
      if (jlat1.eq.ny .and. xlat.lt.xlatn) jlat1 = ny - 1   ! allow extrapolation to xlatn
      jlat2 = jlat1 + 1
      if (jlat1.lt.1 .or. jlat2.gt.ny) then
         isdata = .false.

*     if the orthoweight file repesents differences return zero
*     instead of saying there is no data
*     ---------------------------------------------------------
         if (idiff.eq.1) isdata = .true.
         return
      endif
      xlon = dlon
      if (xlon.lt.lonmin) xlon = xlon + 360.d0 - eps   ! eps to avoid problem when xlon=lonmin
      ilon1 = int((xlon - lonmin)/dx) + 1
      if (ilon1.gt.nx) ilon1 = ilon1 - 1               ! avoid problem for dlon=360.
      ilon2 = ilon1 + 1
      if (ilon2.gt.nx) ilon2 = 1                       ! longitude grid must span globe
      if (ilon1.lt.1 .or. ilon1.gt.nx) then
         write (6,*) 'error in tptide: incorrect longitude input'
         write (6,*) 'mjd,dlon,xlon,ilon1,ilon2',
     .                time/86400.d0,dlon,xlon,ilon1,ilon2
         stop 301
      endif

      if (uv(1,1,1,ilon1,jlat1).eq.undef .and.
     *    uv(1,1,1,ilon2,jlat1).eq.undef .and.
     *    uv(1,1,1,ilon1,jlat2).eq.undef .and.
     *    uv(1,1,1,ilon2,jlat2).eq.undef)  then
         isdata = .false.
         if (idiff.eq.1) isdata = .true.
         return
      endif

*  distances from corners in units of cell size
*  --------------------------------------------
      wx1 = ( dx - (xlon - real(ilon1-1)*dx - lonmin) )/dx
      wx2 = 1.d0 - wx1
      wy1 = ( dy - (xlat - real(jlat1-1)*dy - latmin) )/dy
      wy2 = 1.d0 - wy1
      if (wy2.lt.0.) then           ! north corner weights zero in southern extrapolation
         wy1 = wy1 + wy2
         wy2 = 0.
      endif
      if (wy1.lt.0.) then           ! south corner weights zero in northern extrapolation
         wy2 = wy1 + wy2
         wy1 = 0.
      endif

*  interpolation weights:
*  w1,w2,w3,w4 are for northwest,northeast,southeast,southwest corners.
*  --------------------------------------------------------------------
      w1 = wx1*wy2
      w2 = wx2*wy2
      w3 = wx2*wy1
      w4 = wx1*wy1

*  initialize sum of weights and orthoweights
*  ------------------------------------------
      w = 0.d0
      do 50 i=1,3
      do 50 j=1,2
      u(i,j) = 0.d0
   50 v(i,j) = 0.d0

*     interpolate the orthoweights
*     ----------------------------
      if (uv(1,1,1,ilon1,jlat1).ne.undef) then
         do 100 i=1,3
            u(i,1) = uv(1,i,1,ilon1,jlat1)*w4
            v(i,1) = uv(2,i,1,ilon1,jlat1)*w4
            u(i,2) = uv(1,i,2,ilon1,jlat1)*w4
            v(i,2) = uv(2,i,2,ilon1,jlat1)*w4
  100    continue
         w = w4
      endif
      if (uv(1,1,1,ilon1,jlat2).ne.undef) then
         do 200 i=1,3
            u(i,1) = uv(1,i,1,ilon1,jlat2)*w1 + u(i,1)
            v(i,1) = uv(2,i,1,ilon1,jlat2)*w1 + v(i,1)
            u(i,2) = uv(1,i,2,ilon1,jlat2)*w1 + u(i,2)
            v(i,2) = uv(2,i,2,ilon1,jlat2)*w1 + v(i,2)
  200    continue
         w = w + w1
      endif
      if (uv(1,1,1,ilon2,jlat2).ne.undef) then
         do 300 i=1,3
            u(i,1) = uv(1,i,1,ilon2,jlat2)*w2 + u(i,1)
            v(i,1) = uv(2,i,1,ilon2,jlat2)*w2 + v(i,1)
            u(i,2) = uv(1,i,2,ilon2,jlat2)*w2 + u(i,2)
            v(i,2) = uv(2,i,2,ilon2,jlat2)*w2 + v(i,2)
  300    continue
         w = w + w2
      endif
      if (uv(1,1,1,ilon2,jlat1).ne.undef) then
         do 400 i=1,3
            u(i,1) = uv(1,i,1,ilon2,jlat1)*w3 + u(i,1)
            v(i,1) = uv(2,i,1,ilon2,jlat1)*w3 + v(i,1)
            u(i,2) = uv(1,i,2,ilon2,jlat1)*w3 + u(i,2)
            v(i,2) = uv(2,i,2,ilon2,jlat1)*w3 + v(i,2)
  400    continue
         w = w + w3
      endif
      if (w.gt.wmin) then
         do 500 i=1,3
         do 500 j=1,2
         u(i,j) = u(i,j)/w
  500    v(i,j) = v(i,j)/w
         call tptid1( u,v,time,pseudo,radiation,tide )
      else
         isdata = .false.
         if (idiff.eq.1) then
            isdata = .true.
            tide = 0.0d0
            do i=1,3
               do j=1,2
                  u(i,j) = 0.d0
                  v(i,j) = 0.d0
               enddo
            enddo
         endif
      endif
      return
      end


      subroutine tptid1( u,v,ts,pseudo,radiation,tide )

*  computes the tide at time ts from the orthoweights u,v.
*  this is a kernel routine, meant to be called only by tptide.

      implicit double precision (a-h,o-z)
      dimension   indx(30,2,4),amp(30,2),freq(30,2),
     *            pha(30,2),ct(30,2),st(30,2),at(3),bt(3),
     *            u00(2),u20(2),u21(2),u40(2),u41(2),v41(2),
     *            p(3,2),q(3,2),u(3,2),v(3,2),
     *            shpn(4),phc(4),dpd(4)
      logical     init,radiation,pseudo
      parameter  (rad_factor=0.97d0, rad_phase=5.9d0)
      save
*
      data phc/290.21d0, 280.12d0, 274.35d0, 343.51d0/,
     *     dpd/13.1763965d0,0.9856473d0,0.1114041d0,0.0529539d0/,
     *     u00/0.0298d0,0.0200d0/, u20/0.1408d0,0.0905d0/,
     *     u21/0.0805d0,0.0638d0/, u40/0.6002d0,0.3476d0/,
     *     u41/0.3025d0,0.1645d0/, v41/0.1517d0,0.0923d0/
      data tc/40431744.d2/, tslast/-9.99e10/
      data init/.true./
      data indx                                             !(  l  ,m,n)
     + /-3,-3,-2,-2,-2,-2,-1,-1,-1,-1,-1, 0, 0, 0, 0,       !( 1:15,1,1)
     +   1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3,       !(16:30,1,1)
     +  -3,-3,-2,-2,-2,-1,-1,-1,-1,-1,-1, 0, 0, 0, 0,       !( 1:15,2,1)
     +   0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3,       !(16:30,2,1)
     +   0, 2, 0, 0, 2, 2, 0, 0, 0, 0, 2, 0, 0, 0, 2,       !( 1:15,1,2)
     +  -3,-2,-2, 0, 0, 0, 1, 2,-2, 0, 0,-2, 0, 0, 0,       !(16:30,1,2)
     +   0, 2, 0, 2, 3,-1, 0, 0, 1, 2, 3,-2,-1, 0, 0,       !( 1:15,2,2)
     +   1,-2, 0, 0, 0,-3,-2,-1, 0, 0, 0, 0,-2, 0, 0,       !(16:30,2,2)
     +   2, 0, 1, 1,-1,-1, 0, 0, 0, 2, 0,-1, 1, 1,-1,       !( 1:15,1,3)
     +   0, 0, 0, 0, 0, 0, 0, 0, 1,-1,-1, 0,-2, 0, 0,       !(16:30,1,3)
     +   3, 1, 2, 0, 0, 1, 1, 1, 1,-1,-1, 2, 0, 0, 0,       !( 1:15,2,3)
     +   0, 1,-1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1,-1,-1,       !(16:30,2,3)
     +   0, 0,-1, 0,-1, 0,-2,-1, 0, 0, 0, 0, 0, 1, 0,       !( 1:15,1,4)
     +   0,-1, 0,-1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1,       !(16:30,1,4)
     +   0, 0, 0, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0,-1, 0,       !( 1:15,2,4)
     +   0, 0, 0, 0, 1, 0, 0, 0,-1, 0, 1, 2, 0, 0, 1/       !(16:30,2,4)
      data amp                                                !(  l  ,m)
     + /0.00663,0.00802,0.00947,0.05019,0.0018 ,0.00954,      !( 1: 6,1)
     +  0.00152,0.04946,0.26218,0.00171,0.00343,0.00741,      !( 7:12,1)
     +  0.02062,0.00414,0.00394,0.00713,0.00137,0.122  ,      !(13:18,1)
     +  0.0073 ,0.36874,0.05002,0.00293,0.00524,0.00395,      !(19:24,1)
     +  0.02061,0.00409,0.00342,0.00169,0.01128,0.00723,      !(25:30,1)
     +  0.0018 ,0.00467,0.01601,0.01932,0.0013 ,0.00102,      !( 1: 6,2)
     +  0.00451,0.121  ,0.00113,0.02298,0.00106,0.0019 ,      !( 7:12,2)
     +  0.00218,0.02358,0.63194,0.00193,0.00466,0.01786,      !(13:18,2)
     +  0.00447,0.00197,0.01718,0.29402,0.00305,0.00102,      !(19:24,2)
     +  0.07994,0.02382,0.00259,0.00086,0.00447,0.00195/      !(25:30,2)
      data pha                                                !(  l  ,m)
     + / 90., 90., 90., 90., 90., 90.,270., 90., 90.,270.,    !( 1:10,1)
     +  270.,270.,270.,270.,270., 13.,270., 90., 90.,270.,    !(11:20,1)
     +  270.,-13.,270.,270.,270.,270.,270.,270.,270.,270.,    !(21:30,1)
     +    0.,  0.,  0.,  0., 77.,103.,180.,  0., 77.,  0.,    !( 1:10,2)
     +   77.,180.,103.,180.,  0., 77.,180.,180.,  0.,  0.,    !(11:20,2)
     +  283.,  0.,262.,180.,  0.,  0.,  0.,  0.,  0.,  0./    !(21:30,2)
*
      if (init) then
      init = .false.
      rad = datan(1.d0)/45.d0
      dpld  = 36.d1 - dpd(1) + dpd(2)
      twopi = datan(1.d0)*8.d0
*
      fdpd  = 0.d0
c     write(6,100)
c 100 format(/'   semidiurnal & diurnal tidal potential terms:'/)
      do 201 m=1,2
         fdpd = fdpd + dpld
         do 101 l=1,30
c  the data formerly read here is now in the data statements above
c           read(9,*) (indx(l,m,n),n=1,4), amp(l,m), pha(l,m)
            theta = fdpd
            do 1 n=1,4
               theta = theta + real(indx(l,m,n))*dpd(n)
    1       continue
            freq(l,m) = theta
            omegat = 2.d0*theta*rad
            ct(l,m) = 2.d0*dcos(omegat)
            st(l,m) = 2.d0*dsin(omegat)
            if (radiation .and. m.eq.2 .and. 
     *           (l.eq.21.or.l.eq.22.or.l.eq.23)) then
               amp(l,m) = amp(l,m)*rad_factor
               pha(l,m) = pha(l,m) - rad_phase
               write(6,503) l,m,(indx(l,m,n),n=1,4),freq(l,m),amp(l,m),
     *                      pha(l,m)
            else
c              write(6,502) l,m,(indx(l,m,n),n=1,4),freq(l,m),amp(l,m),
c    *                      pha(l,m)
            endif
  101    continue
c        write(6,*) ' '
  201 continue
  502 format(1x,i3,i4,2x,4i2,f12.6,2(f10.5,f10.1,5x))
  503 format(1x,i3,i4,2x,4i2,f12.6,f10.5,f10.1,3x,'incl.rad.')
      endif

      if (ts.ne.tslast) then
         td = (ts - tc)/864.d2
*    compute 4 principal mean longitudes for a given time td
         do 107 n=1,4
            ph = phc(n) + td*dpd(n)
            shpn(n) = dmod(ph,36.d1)
  107    continue
         fd = td - dint(td)
         e  = 36.d1*fd - shpn(1) + shpn(2)
      endif
*
      otide = 0.d0
      do 300 m=1,2
         if (ts.ne.tslast) then
*
*  compute orthotides p(3),q(3) for given time ts in terms
*  of potential amplitudes at(3), bt(3)
            do 207 k=1,3
               at(k) = 0.d0
               bt(k) = 0.d0
  207       continue
            do 407 l=1,30
               ph = m*e + pha(l,m)
               do 307 n=1,4
                  i = indx(l,m,n)
                  if (i.eq.0) go to 307
                  ph = ph + real(i)*shpn(n)
  307          continue
               theta = ph*rad
               a = amp(l,m)*100.d0
               cr = a*dcos(theta)
               sr = a*dsin(theta)
               crct = cr*ct(l,m)
               crst = cr*st(l,m)
               srct = sr*ct(l,m)
               srst = sr*st(l,m)
               at(1) = at(1) + cr
               bt(1) = bt(1) - sr
               at(2) = at(2) + crct
               bt(2) = bt(2) - srct
               at(3) = at(3) + crst
               bt(3) = bt(3) - srst
  407       continue
            p(1,m) = u00(m)*at(1)
            q(1,m) = u00(m)*bt(1)
            p(2,m) = u20(m)*at(1) - u21(m)*at(2)
            q(2,m) = u20(m)*bt(1) - u21(m)*bt(2)
            if (pseudo) then
*  cartwright and ray pseudo-orthoweights
               p(3,m) = u40(m)*at(1) - u41(m)*at(2) - v41(m)*at(3)
               q(3,m) = u40(m)*bt(1) - u41(m)*bt(2) - v41(m)*bt(3)
            else
*  correct orthoweight definition
               p(3,m) = u40(m)*at(1) - u41(m)*at(2) + v41(m)*at(3)
               q(3,m) = u40(m)*bt(1) - u41(m)*bt(2) + v41(m)*bt(3)
            endif
         endif
         do 30 k=1,3
            otide = otide + u(k,m)*p(k,m) + v(k,m)*q(k,m)
   30    continue
  300 continue
      tide = otide
      tslast = ts
      return
      end


      subroutine admit2 (pseudo,radiation,u,v,x,y,h1,h2,amp,pha)
c
c...to compute admittance and in-phase and quadrature components
c   from orthoweights
c
c...coded by Richard Eanes, UT/CSR, 1993 and 1994
c
c Inputs:
c  pseudo   -- Logical variable, true for CR91 style pseudo-orthoweights
c              Use output from tptide for this
c  radiation - Logical variable, true if the orthoweights are corrected
c              for the radiation potential at T2,S2,R2  (R2 not in this list)
c  u,v       - The ortwhoweights for the diurnal and semidiurnal bands
c
c Outputs:
c  x,y       - The tidal admittances at a list of frequencies spanning
c              the diurnal and semidiurnal bands (see data below)
c  h1,h2     - In-phase and quadrature amplitudes at the same frequencies (cm)
c  amp,pha   - Amplitude and phase at the same list of frequencies. (cm,deg)
c
      implicit double precision (a-h,o-z)
      parameter (nf=30)
      logical first,pseudo,radiation
      character*11 tname
      dimension u(3,2), v(3,2)
      dimension amp(100),pha(100),x(100),y(100),h1(100),h2(100)
      dimension tname(nf),freq(nf),per(nf),hsm(nf),chi(nf)
      dimension i1(0:2),i2(0:2),xp1(nf),xp2(nf),fxy(nf),sp(6,2)
      save

c...the orthotide weight factors (signs with + were wrong in C&R 90 Table A1)
      data ((sp(i,m),i=1,6),m=1,2) /
     . 0.0298,  0.1408, +0.0805,  0.6002, +0.3025,  0.1517,
     . 0.0200,  0.0905, +0.0638,  0.3476, +0.1645,  0.0923/
c
      data (tname(i),freq(i),per(i),hsm(i),chi(i),i=1,nf) /
c              name        f(cycle/d)       period      hs (m)    phase
c long period
     .     'LP   55.565',  0.000147094,    6798.382,    0.0279,   180.0, ! 1
     .     'Sa   56.554',  0.002737779,     365.260,   -0.0049,     0.0, ! 2
     .     'Ssa  57.555',  0.005475819,     182.621,   -0.0310,     0.0, ! 3
     .     'TERa 58.554',  0.008213597,     121.749,   -0.0018,     0.0, ! 4
     .     'Mm   65.455',  0.036291647,      27.555,   -0.0352,     0.0, ! 5
     .     'Mf   75.555',  0.073202203,      13.661,   -0.0666,     0.0, ! 6
     .     'TERm 85.455',  0.109493850,       9.133,   -0.0128,     0.0, ! 7
     .     '93a  93.555',  0.140928587,       7.096,   -0.0020,     0.0, ! 8
c diurnal
     .     '2Q1 125.755',  0.856952384,      -6.904,   -0.0066,   270.0, ! 9
     .     'Q1  135.655',  0.893244048,      -9.213,   -0.0502,   270.0, !10
     .     'O1  145.555',  0.929535695,     -13.841,   -0.2622,   270.0, !11
     .     'M1  155.655',  0.966446233,     -28.297,    0.0206,    90.0, !12
     .     'P1  163.555',  0.997262079,    -221.060,   -0.1220,   270.0, !13
     .     'S1  164.556',  1.000000000,    -559.897,    0.0029,    90.0, !14
     .     'K1  165.555',  1.002737898,    1050.246,    0.3688,    90.0, !15
     .     'PHI1167.555',  1.008213699,     155.580,    0.0053,    90.0, !16
     .     'J1  175.455',  1.039029527,      26.850,    0.0206,    90.0, !17
     .     'OO1 185.555',  1.075940083,      13.485,    0.0113,    90.0, !18
     .     'NU1 195.455',  1.112231730,       9.054,    0.0022,    90.0, !19
c semi-diurnal
     .     '227 227.655',  1.828255527,      -5.704,    0.0047,     0.0, !20
     .     '2N2 235.755',  1.859690264,      -6.950,    0.0160,     0.0, !21
     .     'MU2 237.555',  1.864547173,      -7.193,    0.0193,     0.0, !22
     .     'N2  245.655',  1.895981946,      -9.295,    0.1210,     0.0, !23
     .     'NU2 247.455',  1.900838820,      -9.734,    0.0230,     0.0, !24
     .     'M2  255.555',  1.932273593,     -14.026,    0.6319,     0.0, !25
     .     'L2  265.455',  1.968565204,     -28.566,   -0.0179,   180.0, !26
     .     'T2  272.556',  1.997262163,    -158.475,    0.0172,     0.0, !27
     .     'S2  273.555',  2.000000000,    -279.994,    0.2940,     0.0, !28
     .     'K2  275.555',  2.005475795,     525.123,    0.0800,     0.0, !29
     .     '285 285.455',  2.041767407,      26.181,    0.0045,     0.0/ !30
      data i1 / 1, 9, 20/,  i2 / 8, 19, 30/
      data dt /2.d0/
      data rad_fac /0.97d0/, rad_phase /5.9d0/   ! Ray-Sanchez-Cartrwight radiation factors
      data first /.true./

c...compute constants depending only on frequency
      first = .true.
      if (first) then
         first = .false.
         twopi = 8.d0*atan(1.d0)
         deg = 360.d0/twopi
         do m = 1,2
           mf = (-1)**m
           do i = i1(m),i2(m)
              theta = twopi*freq(i)*dt
              cw = 2.*cos(theta)
              sw = 2.*sin(theta)
              xp1(i) = sp(2,m) - sp(3,m)*cw
              if (pseudo) then
* Cartwright&Ray pseudo-orthoweights
                 xp2(i) = sp(4,m) - sp(5,m)*cw - sp(6,m)*sw
              else
* True orthoweights
                 xp2(i) = sp(4,m) - sp(5,m)*cw + sp(6,m)*sw
              endif
              fxy(i) = mf*abs(hsm(i))*1000.d0   ! scale to return values in mm
c             write (12,200) tname(i),fxy(i)*sp(1,m),fxy(i)*xp1(i),
c    .                                               fxy(i)*xp2(i)
           enddo
         enddo
      endif

      do m = 1,2
         do i = i1(m),i2(m)
            x(i) = sp(1,m)*u(1,m) + xp1(i)*u(2,m) + xp2(i)*u(3,m)
            y(i) = sp(1,m)*v(1,m) + xp1(i)*v(2,m) + xp2(i)*v(3,m)
            h1(i) =  fxy(i)*x(i)
            h2(i) = -fxy(i)*y(i)
            amp(i) = sqrt(h1(i)**2 + h2(i)**2)
            if (amp(i).ne.0.) then
               pha(i) = atan2(h2(i),h1(i))*deg
               if (pha(i).lt.0.) pha(i) = pha(i) + 360.
            else
               pha(i) = 0.
            endif
c...correct T2 and S2 for radiation potential
            if (radiation .and. (i.eq.27 .or. i.eq.28) ) then
               amp(i) = amp(i)*rad_fac
               pha(i) = pha(i) + rad_phase
               h1(i) = amp(i)*cos(pha(i)/deg)
               h2(i) = amp(i)*sin(pha(i)/deg)
            endif
c           write (6,100) i,m,tname(i),hsm(i),freq(i),x(i),y(i),
c    .                    h1(i),h2(i),amp(i),pha(i)
        enddo
      enddo

  100 format (1x,2i2,1x,a11,1x,f8.4,f8.5,1x,4f9.4,1x,2f8.2)
  200 format (a11,3x,3f15.6)

      return
      end


      subroutine lpeqmt( ts, dlat, tlp )
*
*  function - computes the long-period equilibrium ocean tides.
*
*  arguments -
*     name      type  i/o               description
*     ----      ----  ---               -----------
*     ts         d     i   modified julian day, in seconds, denoting
*                          time at which tide is to be computed.
*
*     dlat       d     i   latitude in degrees (positive north) for the
*                          position at which tide is computed.
*
*     tlp        d     o   computed long-period tide, in centimeters.
*
*
*  processing logic -
*     fifteen tidal spectral lines from the cartwright-tayler-edden
*     tables are summed over to compute the long-period tide.
*
*  technical references -
*     cartwright & tayler, geophys. j. r.a.s., 23, 45, 1971.
*     cartwright & edden, geophys. j. r.a.s., 33, 253, 1973.
*
*  history -
*    version   date       programmer         change description
*    -------   ----       ----------         ------------------
*      1.0   11/27/90     d. cartwright    documented by r. ray
*
*
      implicit double precision (a-h,o-z)
      dimension  phc(4),dpd(4),shpn(4)
      save
*
      parameter (rad= 0.017453292519943 d0)
      parameter (psol=283.d0*rad)       ! solar perigee held constant
      data phc/290.21d0,280.12d0,274.35d0,343.51d0/,
     *     dpd/13.1763965d0,0.9856473d0,0.1114041d0,0.0529539d0/,
     *     tc/40431744.d2/
*
*     compute 4 principal mean longitudes in radians at time td
*     ---------------------------------------------------------
      td = (ts - tc)/864.d2             ! days past 1 jan 87 00:00
      do 1 n=1,4
         ph = phc(n) + td*dpd(n)
         shpn(n) = rad*dmod(ph,36.d1)
    1 continue
*
*     assemble long-period tide potential from 15 cte terms > 1 mm.
*     nodal term is included but not the constant term.
*     -------------------------------------------------
      zlp = 2.79d0 * dcos( shpn(4) )
     *     -0.49d0 * dcos( shpn(2) - psol )
     *     -3.10d0 * dcos( 2.d0*shpn(2) )
      ph = shpn(1)
      zlp = zlp - 0.67d0 * dcos( ph - 2.d0*shpn(2) + shpn(3) )
     *          -(3.52d0 - 0.46d0*dcos(shpn(4))) * dcos( ph - shpn(3) )
      ph = ph + shpn(1)
      zlp = zlp - 6.66d0 * dcos( ph )
     *          - 2.76d0 * dcos( ph + shpn(4) )
     *          - 0.26d0 * dcos( ph + 2.d0*shpn(4) )
     *          - 0.58d0 * dcos( ph - 2.d0*shpn(2) )
     *          - 0.29d0 * dcos( ph - 2.d0*shpn(3) )
      ph = ph + shpn(1)
      zlp = zlp - 1.27d0 * dcos( ph - shpn(3) )
     *          - 0.53d0 * dcos( ph - shpn(3) + shpn(4) )
     *          - 0.24d0 * dcos( ph - 2.d0*shpn(2) + shpn(3) )
*
*     multiply by gamma_2 * sqrt(5/4 pi) * p20(lat)
*     ---------------------------------------------
      s = dsin(dlat*rad)
      tlp = 0.437d0*zlp*(1.5d0*s*s - 0.5d0)
      return
      end


      subroutine kalday (utc,mdyhms)
c
      implicit double precision (a-h,o-z)
      dimension mdyhms(6),iptime(5)
c
c...compute the hour, minute, and second.
c
      t = utc
    5 jd = t
      xtime = mod((t-jd)*24.0d0+12.0d0,24.0d0)
      iptime(4) = xtime
      xtime = (xtime-iptime(4))*60.0d0
      iptime(5) = xtime
      ptimes = (xtime-iptime(5))*60.0d0
c
c...compute the year, month, and day.
c
      if (iptime(4).lt.12) jd = jd+1
      lx = jd+68569
      nx = 4*lx/146097
      lx = lx-(146097*nx+3)/4
      iptime(3) = 4000*(lx+1)/1461001
      lx = lx-1461*iptime(3)/4+31
      iptime(1) = 80*lx/2447
      iptime(2) = lx-2447*iptime(1)/80
      lx = iptime(1)/11
      iptime(1) = iptime(1)+2-12*lx
      iptime(3) = 100*(nx-49)+iptime(3)+lx
      do 10 i = 1,5
         mdyhms(i) = iptime(i)
   10 continue
      mdyhms(3) = mdyhms(3)-1900
      mdyhms(6) = ptimes+0.5d0
c
      if (mdyhms(4).eq.23 .and. mdyhms(5).eq.59 .and.
     .  mdyhms(6).eq.60) then
        t = t + 1.0d-6
        go to 5
      endif
c
      return
c
      end
      subroutine julday (itimes,stime,utc)
c
c...purpose:  to convert between calendar day and julian date (utc).
c
c   output:
c       utc      double precision julian date
c
c   input:
c      itimes    integer array containing month, day, year, hour, minute
c      stime     floating point seconds
c
      implicit double precision (a-h,o-z)
      dimension itimes(5)
c
      jd = itimes(2)-32075+1461*(itimes(3)+4800+(itimes(1)-14)/12)/4+367
     *   *(itimes(1)-2-(itimes(1)-14)/12*12)/12-3*((itimes(3)+4900+
     *   (itimes(1)-14)/12)/100)/4
c
      utc = dble(float(jd))-0.5d+0+dble(float(itimes(4)))/24.0d+0+
     *   dble(float(itimes(5)))/1440.0d+0+dble(stime)/86400.0d+0
c
      return
c
c.....end julday
c
      end
