cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      SUBROUTINE DWOTIDE(LAT,LON,TIME,ITYPE,TIDE1,TIDE2,ISDATA)
c.....This subroutien serves to read the orthoweight tables that define
c.....the Desai and Wahr (1995) ocean tide model derived from
c.....TOPEX/POSEIDON and extended globally with Schwiderski's
c.....ocean tide model
c..... REFERENCES:
c......  Groves and Reynolds [1975], J. Geophys. Res. v80, 4131-4138
c.....   Cartwright and Taylor [1971], J. Geophys. Res. v23. 45-74
c.....   Desai and Wahr [1995], J. Geophys. Res., in press
c.....The routine uses bilinear interpolation to determine the
c.....orthoweights at the specified coordinates.
c.....
c.....ADDITIONAL SUBROUTINES CALLED FROM THIS SUBROUTINE
c.....  ORTHOWEIGHTRES - Defines orthotide constants and lag interval
c.....  READAMP - Reads tidal components with Doodson numbers
c.....              tide potential amplitudes and computes
c.....              frequencies and astronomical phase angles
c.....  DWOUT   - Computes one of 4 Options for the ocean tides
c.....              including total time dependent ocean tide heights,
c.....              admittance, in phase and quadrature components,
c.....              and amplitude and Greenwich phase
c.....
c.....INPUT DATA FILES (Should be opened before calling this routine)
c.....  UNIT=50 - tidepot_dw.amp (with specified tidal components to include)
c.....  UNIT=60 - Orthoweights files 
c.....        eg. cuelastic_dw3.2.90 for elastic ocean tide (ocean+load)
c.....        eg. cuocean_dw3.2.90 for pure ocean tide 
c.....        eg. cuload_dw3.2.90 for pure load tide 
c.....
c.....SUBROUTINE INPUTS:
c.....  LAT  - North Latitude (degrees)
c.....  LON  - East Longitude (degrees)
c.....  TIME - Time in MJD
c.....  ITYPE- Specifies one of four options for output
c.....    ITYPE = 1 => Total time dependent ocean tide
c.....    ITYPE = 2 => Real and imaginary admittance of 1 component
c.....    ITYPE = 3 => In phase and Quadrature components of 1 component
c.....    ITYPE = 4 => Amplitude and Greenwich Phase of 1 component
c.....  NOTE: if ITYPE=2,3,4 ONLY ONE tidal component can be specified
c.....        in tidepot_dw.amp
c.....
c.....SUBROUTINE OUTPUT:
c.....  TIDE1 - Returns
c.....    ITYPE = 1, Total time dependent ocean tide height (cm)
c.....    ITYPE = 2, Real Admittance (%)
c.....    ITYPE = 3, In Phase Component (cm)
c.....    ITYPE = 4, Amplitude (cm)
c.....  TIDE2 - Returns
c.....    ITYPE = 1, Zero (Empty field)
c.....    ITYPE = 2, Imaginary Admittance (%)
c.....    ITYPE = 3, Quadrature Component (cm)
c.....    ITYPE = 4, Greenwich Phase (degrees)
c.....  ISDATA - Logical variable returns false if not valid
c.....           model at specified (LAT,LON)
c.....
c..... AUTHOR: SHAILEN D. DESAI
c..... DATE  : 09/13/95
c.....
      IMPLICIT NONE
      DOUBLE PRECISION LAT, LON, TIME, TIDE1, TIDE2
      INTEGER IFCN, ITYPE
      INTEGER NTIDE, CMAX
      INTEGER NORTHO, IMAX, JMAX
      INTEGER NPAR2, NPAR1
      PARAMETER (IFCN=1)
      PARAMETER (NTIDE=7, CMAX=162)
      PARAMETER (NORTHO=22, IMAX=360, JMAX=180)
      PARAMETER (NPAR2=2)
      PARAMETER (NPAR1=(NORTHO-5*NPAR2)/2)
      DOUBLE PRECISION  U(7,5,0:2), V(7,5,0:2), TAU
      DOUBLE PRECISION AMP(NTIDE,CMAX), FRE(NTIDE,CMAX)
      DOUBLE PRECISION UNDEF
      DOUBLE PRECISION LATMIN, LATMAX, LONMIN, LONMAX
      DOUBLE PRECISION ALPHA(IMAX,JMAX,NORTHO), TEMP(NORTHO)
      DOUBLE PRECISION LAT, LON
      DOUBLE PRECISION DX, DY
      DOUBLE PRECISION LON1, LAT1, WLON1, WLON2, WLAT1, WLAT2
      DOUBLE PRECISION W11, W12, W21, W22
      DOUBLE PRECISION W, ALPIJ(NTIDE,NPAR1)
      INTEGER D(NTIDE,CMAX,6), NMAX(NTIDE), NTOT
      INTEGER IB, JB
      INTEGER IJ, K, N
      INTEGER NX, NY
      INTEGER ILON1, ILON2, ILAT1, ILAT2
      INTEGER KMAX, K0, KK
      CHARACTER*80 TITLE1, TITLE2
      LOGICAL ISDATA
      LOGICAL INIT
      DATA INIT / .TRUE. /
      DATA UNDEF / 99999.0D0 /

c.....First check ITYPE within bounds
      IF ((ITYPE .LT. 1) .OR. (ITYPE .GT. 4)) THEN
	WRITE(6,*) '***FATAL ERROR***'
	WRITE(6,*) 'ITYPE MUST BE WITHIN 1 AND 4'
	WRITE(6,*) '***FATAL ERROR***'
	STOP
      ENDIF
c.....
      IF (INIT) THEN
	INIT=.FALSE.
c.....
c.....Read orthotide constants
	CALL ORTHOWEIGHTRES(U,V,TAU)
c.....
c.....Read Specified tidal components
        CALL READAMP(NTIDE,CMAX,IFCN,AMP,D,FRE,NMAX,NTOT)
	WRITE(6,*)
	WRITE(6,*) 'Total number of specified tidal components =',NTOT
	WRITE(6,*)
	IF ((ITYPE .NE. 1) .AND. (NTOT .NE. 1))THEN
	  WRITE(6,*) '*** FATAL ERROR!!***'
	  WRITE(6,*) 'MUST SPECIFY ONLY 1 TIDAL COMPONENT IF ITYPE>1'
	  WRITE(6,*) '*** FATAL ERROR!!***'
	  STOP
        ENDIF
c.....
c.....Read orthoweights from tables
        READ(60,1400) TITLE1, TITLE2
	READ(60,1500) NX,NY,LATMIN,LATMAX,LONMIN,LONMAX
	IF ((NX .GT. IMAX).OR. (NY .GT. JMAX)) THEN
	  WRITE(6,*) 'INCORRECT DIMENSIONS ON IMAX AND JMAX'
	  WRITE(6,*) NX,NY,IMAX,JMAX
	  STOP
        ENDIF
	WRITE(6,1400) TITLE1,TITLE2
	WRITE(6,1500) NX,NY,LATMIN,LATMAX,LONMIN,LONMAX
	DX=(LONMAX-LONMIN)/DBLE(NX-1)
	DY=(LATMAX-LATMIN)/DBLE(NY-1)
	WRITE(6,*) 'ORTHOWEIGHT TABLES IN INTERVALS OF ',DX,DY
c.....Initialize all orthoweights to UNDEF
	DO 10 IB=1,IMAX
	  DO 20 JB=1,JMAX
	    DO 30 K=1,NORTHO
	      ALPHA(IB,JB,K)=UNDEF
30          CONTINUE
20        CONTINUE
10      CONTINUE
c.....
	DO 40 IJ=1,NX*NY
	  READ(60,1200,END=51) JB,IB,LAT1,LON1,(TEMP(K),K=1,NORTHO)
	  DO 50 K=1,NORTHO
	    ALPHA(IB,181-JB,K)=TEMP(K)
50        CONTINUE
40      CONTINUE
51      CONTINUE
        CLOSE(60)
c.....
1200  FORMAT(2I4,2F6.1,22F10.5)
1400  FORMAT(A80,/,A80)
1500  FORMAT(2I4,4F8.3)
      ENDIF
c.....
c.....Find nearest lower left corner of grid points surrounding
c.....specified point (LON,LAT)
c.....
c.....  (LON1,LAT1+DY)                (LON1+DX,LAT1+DY)
c.....       X--------------------------------X
c.....       |                                |
c.....       |                                |
c.....       |        +(LON,LAT)              |
c.....       |                                |
c.....       X--------------------------------X
c.....  (LON1,LAT1)                   (LON1+DX,LAT1)
c.....
      IF (LON .LT. LONMIN) THEN
	LON=LON+360.0D0
      ENDIF
      ILON1=INT((LON-LONMIN)/DX)+1
      ILON2=ILON1+1
      IF (ILON2 .GT. NX) THEN
	ILON2=ILON2-NX
      ENDIF
      IF ((ILON1 .LT. 1) .OR. (ILON1 .GT. NX)) THEN
        WRITE(6,*) 'ERROR WITH LON. GRID VALUES',ILON1,LON
	STOP
      ENDIF
      ILAT1=INT((LAT-LATMIN)/DY)+1
      ILAT2=ILAT1+1
      IF ((ILAT1 .LT. 1) .OR. (ILAT1 .GT. NY)) THEN
	WRITE(6,*) 'ERROR WITH LAT. GRID VALUES',ILAT1,LAT
	STOP
      ENDIF
      LON1=DBLE(ILON1-1)*DX+LONMIN
      LAT1=DBLE(ILAT1-1)*DY+LATMIN
c.....Compute weight for each corner
      WLON2=(LON-LON1)/DX
      WLON1=1.0D0-WLON2
      WLAT2=(LAT-LAT1)/DY
      WLAT1=1.0D0-WLAT2
      W11=WLON1*WLAT1
      W12=WLON1*WLAT2
      W21=WLON2*WLAT1
      W22=WLON2*WLAT2
c.....
c.....Compute the interpolated orthoweights
      DO 60 N=1,NTIDE
	IF (N .LE. 2) THEN
	  KMAX=NPAR1
	  K0=(N-1)*NPAR1
        ELSE
	  KMAX=NPAR2
	  K0=2*NPAR1+(N-3)*NPAR2
        ENDIF
	DO 70 K=1,KMAX
	  KK=K0+K
	  IF (ALPHA(ILON1,ILAT1,1) .NE. UNDEF) THEN
	    ALPIJ(N,K)=ALPHA(ILON1,ILAT1,KK)*W11
          ELSE
	    ALPIJ(N,K)=0.0D0
          ENDIF
	  IF (ALPHA(ILON1,ILAT2,1) .NE. UNDEF) THEN
	    ALPIJ(N,K)=ALPIJ(N,K)+ALPHA(ILON1,ILAT2,KK)*W12
          ENDIF
	  IF (ALPHA(ILON2,ILAT1,1) .NE. UNDEF) THEN
	    ALPIJ(N,K)=ALPIJ(N,K)+ALPHA(ILON2,ILAT1,KK)*W21
          ENDIF
	  IF (ALPHA(ILON2,ILAT2,1) .NE. UNDEF) THEN
	    ALPIJ(N,K)=ALPIJ(N,K)+ALPHA(ILON2,ILAT2,KK)*W22
          ENDIF
70      CONTINUE
60    CONTINUE
      IF (ALPHA(ILON1,ILAT1,1) .NE. UNDEF) THEN
	W=W11
      ELSE
	W=0.0D0
      ENDIF
      IF (ALPHA(ILON1,ILAT2,1) .NE. UNDEF) THEN
	W=W+W12
      ENDIF
      IF (ALPHA(ILON2,ILAT1,1) .NE. UNDEF) THEN
	W=W+W21
      ENDIF
      IF (ALPHA(ILON2,ILAT2,1) .NE. UNDEF) THEN
	W=W+W22
      ENDIF
c.....
c.....Only compute interpolated orthoweight if W > 0.3
      IF (W .GT. 0.3D0) THEN
	ISDATA=.TRUE.
        DO 80 N=1,NTIDE
	  IF (N .LE. 2) THEN
	    KMAX=NPAR1
          ELSE
	    KMAX=NPAR2
          ENDIF
	  DO 90 K=1,KMAX
	    ALPIJ(N,K)=ALPIJ(N,K)/W
90        CONTINUE
80      CONTINUE
      ELSE
	ISDATA=.FALSE.
	TIDE1=0.0D0
	TIDE2=0.0D0
      ENDIF
c.....
c.....Compute specified output of tide
      IF (ISDATA) THEN
        CALL DWOUT(U,V,TAU,AMP,D,FRE,NMAX,NTIDE,CMAX,ALPIJ,NPAR1,
     &NPAR2,ITYPE,TIME,TIDE1,TIDE2)
      ENDIF



      RETURN
      END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c.....This subroutine reads the specified tidal components
c.....from tidepot_dw.amp and saves their respctive Doodson
c.....numbers, amplitudes and frequencies
c.....
c..... ADDITIONAL SUBROUTINES called from here
c.....   FCNRES - Computes FCN resonance function
c..... INPUT -
c.....   NTIDE- Number of tidal bands
c.....   CMAX - Maximum number of components in tidal bands
c.....   IFCN - IFCN=0 Ignore FCN resonance
c.....          IFCN=1 Include FCN resonance (So redefine amplitudes)
c..... OUTPUT-
c.....   AMP  - Tidal amplitudes (After multiplying by resonance function
c.....          if IFCN=1
c.....   D    - Doodson numbers of specified tidal components
c.....   FRE  - Frequency of tidal components
c.....   NMAX - Saves the number of tidal components in each tidal band
c.....   NTOT - Checks total number of components from all bands
c.....
      SUBROUTINE READAMP(NTIDE,CMAX,IFCN,AMP,D,FRE,NMAX,NTOT)
      INTEGER NTIDE, CMAX, IFCN, D(NTIDE,CMAX,6), NMAX(NTIDE)
      DOUBLE PRECISION AMP(NTIDE,CMAX), FRE(NTIDE,CMAX)
      DOUBLE PRECISION HCTE, RES
      DOUBLE PRECISION OME0(5)
      INTEGER IFLAG, DT(6)
      INTEGER L, M, NT, K
      CHARACTER*20 BAND
      DATA OME0 / 13.17639648D0, 0.98564736D0, 0.11140353D0,
     &             0.05295377D0, 0.00004706D0/

      NTOT=0
      DO 10 N=1,NTIDE
	READ(50,1000) L, M, NT, BAND
	NMAX(N)=0
	DO 20 I=1,NT
	  READ(50,1100) IFLAG,(DT(K),K=1,6),HCTE
	  IF (IFLAG .EQ. 1) THEN
	    NMAX(N)=NMAX(N)+1
	    DO 30 K=1,6
	      D(N,NMAX(N),K)=DT(K)
30          CONTINUE
c.....Convert amplitude to cm
	    AMP(N,NMAX(N))=HCTE*100.0D0
c.....Compute tidal frequencies
	    FRE(N,NMAX(N))=DT(1)*(360.0D0-OME0(1)+OME0(2))/360.0D0
	    DO 40 K=2,6
	      FRE(N,NMAX(N))=FRE(N,NMAX(N))+
     &DBLE(DT(K))*OME0(K-1)/360.0D0
40          CONTINUE
c.....Compute FCN resonance function if necessary
	    IF (IFCN .EQ. 1) THEN
	      IF ((L .EQ. 2) .AND. (M .EQ. 1)) THEN
                CALL FCNRES(FRE(N,NMAX(N)),RES)
		AMP(N,NMAX(N))=AMP(N,NMAX(N))*RES
              ENDIF
            ENDIF
          ENDIF
20      CONTINUE
	WRITE(6,*) 'Including ',NMAX(N),' Terms from ',BAND
	NTOT=NTOT+NMAX(N)
10    CONTINUE
      CLOSE(50)

1000  FORMAT(I1,1X,I1,1X,I3,A20)
1100  FORMAT(I1,1X,6I2,1X,F8.5)
      RETURN
      END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c.....This function computes the FCN resonance function
c.....
      SUBROUTINE FCNRES(FREQ,RES)
      IMPLICIT NONE
      DOUBLE PRECISION FREQ, RES
      DOUBLE PRECISION H210, H21O1, DH21
      DOUBLE PRECISION K210, K21O1, DK21
      DOUBLE PRECISION FREO1, FRESD, FREFCN
      DOUBLE PRECISION GAMMA0
      PARAMETER (H210=0.609D0,     K210=0.302D0)
      PARAMETER (H21O1=0.603D0,    K21O1=0.298D0)
      PARAMETER (DH21=-0.00246D0,  DK21=-0.00123D0)
      PARAMETER (FREO1=0.9295357067D0)
      PARAMETER (FRESD=1.0027379093D0)
      PARAMETER (FREFCN=-(-1.0D0-1/433.2D0)*FRESD)
      PARAMETER (GAMMA0=1.0D0+K210-H210)
      DOUBLE PRECISION H21, K21, GAMMA

c.....Body Love numbers
      H21=H21O1+DH21*(FREQ-FREO1)/(FREFCN-FREQ)
      K21=K21O1+DK21*(FREQ-FREO1)/(FREFCN-FREQ)
      GAMMA=1.0D0+K21-H21
c.....Resonance function
      RES=GAMMA/GAMMA0

      RETURN
      END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c.....Subroutine computes tidal output in form specified by
c.....user 
c..... ADDITIONAL SUBROUTINES called from here
c.....  FINDPHA - COmputes astronomical phase angles
c.....  INPUT :
c.....    U,V - Orthotide constants
c.....    TAU - Lag interval of orthotides
c.....    AMP - Tide potential amplitudes in cm (Includes resonance if used)
c.....    D   - Doodson numbers of tidal components
c.....    FRE - Astronomical tidal frequencies
c.....   NMAX - Number of components specified in each tidal band
c.....   NTIDE- Number of tidal bands
c.....   CMAX - Maximum possible components in each band
c.....   AMPIJ- Interpolated value of orthoweights
c.....   NPAR1- Number of orthotides in diurnal and semidiurnal bands
c.....   NPAR2- Number of orthotides in each of remaining bands
c.....
c.....   TIME - Time in MJD (Only needed for ITYPE=1
c.....   ITYPE- Specifies one of four options for output
c.....    ITYPE = 1 => Total time dependent ocean tide
c.....    ITYPE = 2 => Real and imaginary admittance of 1 component
c.....    ITYPE = 3 => In phase and Quadrature components of 1 component
c.....    ITYPE = 4 => Amplitude and Greenwich Phase of 1 component
c.....  NOTE: if ITYPE=2,3,4 ONLY ONE tidal component can be specified
c.....        in tidepot_dw.amp
c.....
c.....SUBROUTINE OUTPUT:
c.....  TIDE1 - Returns
c.....    ITYPE = 1, Total time dependent ocean tide height (cm)
c.....    ITYPE = 2, Real Admittance (%)
c.....    ITYPE = 3, In Phase Component (cm)
c.....    ITYPE = 4, Amplitude (cm)
c.....  TIDE2 - Returns
c.....    ITYPE = 1, Zero (Empty field)
c.....    ITYPE = 2, Imaginary Admittance (%)
c.....    ITYPE = 3, Quadrature Component (cm)
c.....    ITYPE = 4, Greenwich Phase (degrees)
c.....
      SUBROUTINE DWOUT(U,V,TAU,AMP,D,FRE,NMAX,NTIDE,CMAX,ALPIJ,
     &NPAR1,NPAR2,ITYPE,TIME,TIDE1,TIDE2)
      IMPLICIT NONE
      INTEGER NTIDE, CMAX, NPAR1, NPAR2, ITYPE
      INTEGER D(NTIDE,CMAX,6), NMAX(NTIDE)
      DOUBLE PRECISION U(7,5,0:2), V(7,5,0:2), TAU
      DOUBLE PRECISION AMP(NTIDE,CMAX), FRE(NTIDE,CMAX)
      DOUBLE PRECISION ALPIJ(NTIDE,NPAR1)
      DOUBLE PRECISION TIDE1, TIDE2, TIME
      DOUBLE PRECISION PI, RAD, PHA, PHI(6)
      DOUBLE PRECISION C1, S1, C2, S2, CT, ST
      DOUBLE PRECISION X, Y, COAMP, PHASE, ARGZ
      DOUBLE PRECISION P(6)
      INTEGER N, NB, KMAX, K, KK

      PI=4.0D0*DATAN(1.0D0)
      RAD=PI/180.0D0

c.....Compute total tide height if required
      IF (ITYPE .EQ. 1) THEN
c.....Compute astronimcal phase angles at given time in MJD
        CALL FINDPHASE(TIME,PHI)
	TIDE1=0.0D0
	TIDE2=0.0D0
	DO 10 N=1,NTIDE
	  IF (N .LE. 2) THEN
	    KMAX=NPAR1
          ELSE
	    KMAX=NPAR2
          ENDIF
	  DO 20 K=1,KMAX
	    P(K)=0.0D0
20        CONTINUE
	  DO 30 K=1,NMAX(N)
c.....Compute tidal phase angles
	    IF (N .EQ. 1) THEN
	      PHA=DBLE(D(N,K,1))*PHI(1)-90.0D0
            ELSE
	      PHA=DBLE(D(N,K,1))*PHI(1)
            ENDIF
	    DO 40 KK=2,6
	      PHA=PHA+DBLE(D(N,K,KK))*PHI(KK)
40          CONTINUE
            PHA=DMOD(PHA,360.0D0)
	    CT=AMP(N,K)*DCOS(PHA*RAD)
	    ST=AMP(N,K)*DSIN(PHA*RAD)
	    C1=2.0D0*DCOS(2.0D0*PI*FRE(N,K)*TAU)
	    S1=2.0D0*DSIN(2.0D0*PI*FRE(N,K)*TAU)
            C2=2.0D0*DCOS(2.0D0*PI*FRE(N,K)*2.0D0*TAU)
            S2=2.0D0*DSIN(2.0D0*PI*FRE(N,K)*2.0D0*TAU)
	    P(1)=P(1)  + U(N,1,0)*CT
	    P(2)=P(2)  + U(N,1,0)*ST
	    IF (KMAX .GT. 2) THEN
	      P(3)=P(3)  +(U(N,3,0)+U(N,3,1)*C1)*CT
	      P(4)=P(4)  +(U(N,3,0)+U(N,3,1)*C1)*ST
	      P(5)=P(5)  +(U(N,5,0)+U(N,5,1)*C1+V(N,5,1)*S1)*CT
	      P(6)=P(6)  +(U(N,5,0)+U(N,5,1)*C1+V(N,5,1)*S1)*ST
            ENDIF
30        CONTINUE
          DO 50 K=1,KMAX
	    TIDE1=TIDE1+ALPIJ(N,K)*P(K)
50        CONTINUE
10      CONTINUE
      ELSE
	DO 100 N=1,NTIDE
	  IF (NMAX(N) .GT. 0) THEN
	    NB=N
          ENDIF
100     CONTINUE
	IF (NB .LE. 2) THEN
	  KMAX=NPAR1
        ELSE
	  KMAX=NPAR2
        ENDIF
	C1=2.0D0*DCOS(2.0D0*PI*FRE(NB,1)*TAU)
	S1=2.0D0*DSIN(2.0D0*PI*FRE(NB,1)*TAU)
        C2=2.0D0*DCOS(2.0D0*PI*FRE(NB,1)*2.0D0*TAU)
        S2=2.0D0*DSIN(2.0D0*PI*FRE(NB,1)*2.0D0*TAU)
	X=ALPIJ(NB,1)*U(NB,1,0)
	Y=-ALPIJ(NB,2)*U(NB,1,0)
	IF (KMAX .GT. 2) THEN
	  X=X+ALPIJ(NB,3)*(U(NB,3,0)+U(NB,3,1)*C1)
	  Y=Y-ALPIJ(NB,4)*(U(NB,3,0)+U(NB,3,1)*C1)
	  X=X+ALPIJ(NB,5)*(U(NB,5,0)+U(NB,5,1)*C1+V(NB,5,1)*S1)
	  Y=Y-ALPIJ(NB,6)*(U(NB,5,0)+U(NB,5,1)*C1+V(NB,5,1)*S1)
        ENDIF
	COAMP=DABS(AMP(NB,1)*DSQRT(X**2.0D0+Y**2.0D0))
	IF (X .EQ. 0.0D0) THEN
          IF (Y .EQ. 0.0D0) THEN
	    ARGZ=0.0D0
          ELSEIF (Y .GT. 0.0D0) THEN
	    ARGZ=PI/2.0D0
          ELSE
	    ARGZ=3.0D0*PI/2.0D0
          ENDIF
        ELSE
	  ARGZ=DATAN2(Y,X)
        ENDIF
	IF (NB .EQ. 2) THEN
	  PHASE=-ARGZ
        ELSE
	  PHASE=PI-ARGZ
        ENDIF
c.....Compute admittance in percent
	IF (ITYPE .EQ. 2) THEN
	  TIDE1=X*100.0D0
	  TIDE2=Y*100.0D0
c.....Compute in phase and quadrature components in cm
        ELSEIF (ITYPE .EQ. 3) THEN
	  TIDE1=COAMP*DCOS(PHASE)
	  TIDE2=COAMP*DSIN(PHASE)
c.....Compute coamplitude and cophase in cm and degrees
        ELSEIF (ITYPE .EQ. 4) THEN
	  TIDE1=COAMP
	  TIDE2=PHASE/RAD
	  IF (TIDE2 .LT. 0.0D0) THEN
	    TIDE2=TIDE2+360.0D0
          ELSEIF (TIDE2 .GT. 360.0D0) THEN
	    TIDE2=TIDE2-360.0D0
          ENDIF
        ENDIF
      ENDIF


      RETURN
      END
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c.....This subroutine computes the astronomical phases
c.....Returns phase in degrees
c.....Astronomical phases use Cartwright and Edden [1971]
c.....convention as specified in TOPEX/POSEIDON Solid Earth
c.....Tide Algorithm specifications
c.....Reference time is May 22.5 1960
c.....Phases computed as cubic functions of time
c..... INPUT:  
c.....   MJD    - Specify time in MJD (days)
c..... OUPUT: 
c.....   PHI(6) - Six Astronomcal phase angles in degrees
c.....
      SUBROUTINE FINDPHASE(MJD,PHI)
      IMPLICIT NONE
      DOUBLE PRECISION MJD0
c.....REFERENCE EPOCH: T0 = May 22.5 1960
      PARAMETER (MJD0=37076.5D0)
      DOUBLE PRECISION MJD, PHI(6)
      DOUBLE PRECISION TD, TY, TFRAC
      DOUBLE PRECISION OME0(5), PHI0(5), CPHA(5), DPHA(5)
      INTEGER K
      DATA OME0 / 13.17639648D0, 0.98564736D0, 0.11140353D0,
     &             0.05295377D0, 0.00004706D0/
      DATA PHI0 / 22.22101D0, 60.11923D0, 271.56503D0, 
     &           188.82048D0,282.25919D0/
      DATA CPHA / -0.001133D0, 0.000303D0, -0.010325D0,
     &            -0.002078D0, 0.000450D0/
      DATA DPHA /  0.0000019D0, 0.0000000D0, -0.0000120D0,
     &             0.0000000D0, 0.0000030D0/

c.....
c.....Time after reference epoch (days)
      TD=MJD-MJD0
c.....Time after reference epoch (centuries)
      TY=TD/36525.0D0
c.....Fraction of solar day at any given time
      TFRAC=MJD-DBLE(INT(MJD))
c.....Compute astronomical phase angles at given times
      DO 10 K=2,6
	PHI(K)=PHI0(K-1)+OME0(K-1)*TD
     &+CPHA(K-1)*TY*TY+DPHA(K-1)*TY*TY*TY
	PHI(K)=DMOD(PHI(K),360.0D0)
10    CONTINUE
      PHI(1)=360.0D0*TFRAC-PHI(2)+PHI(3)
      PHI(1)=DMOD(PHI(1),360.0D0)

      RETURN
      END
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c.....Subroutine provides orthotide constants with FCN resonance
c.....OUTPUT :
c.....   U,V - Orthotide constants
c.....   TAU - Lag interval for orthotides
c.....
      SUBROUTINE ORTHOWEIGHTRES(U,V,TAU)
      IMPLICIT NONE
      DOUBLE PRECISION U(7,5,0:2), V(7,5,0:2), TAU
      INTEGER I, J, K

c.....
c.....Lag interval
      TAU=2.0D0
c.....
c.....Orthotide constants
      DO 10 I=1,7
	DO 20 J=1,5
	  DO 30 K=0,2
	    U(I,J,K)=0.0d0
	    V(I,J,K)=0.0d0
30        CONTINUE
20      CONTINUE
10    CONTINUE
      U(1,1,0)= 0.028486D0
      U(1,3,0)= 0.138607D0
      U(1,3,1)=-0.078478D0
      U(1,5,0)= 0.594826D0
      U(1,5,1)=-0.300076D0
      V(1,5,1)= 0.150138D0
      U(2,1,0)= 0.019805D0
      U(2,3,0)= 0.089759D0
      U(2,3,1)=-0.063291D0
      U(2,5,0)= 0.344734D0
      U(2,5,1)=-0.163116D0
      V(2,5,1)= 0.091581D0
      U(3,1,0)= 0.392441D0
      U(4,1,0)= 0.195191D0
      U(5,1,0)= 0.982441D0
      U(6,1,0)= 2.874418D0
      U(7,1,0)= 0.456492D0

      RETURN
      END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
