       SUBROUTINE phase(alon,alat,ph1,ph2)
C     DOUBLE PRECISION VERSION FOR USE ON SUN
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (IM=1801,JM=729)
      PARAMETER (NT=3,NS=2,NC=30)
      PARAMETER (SMALL = 1.0D-8)
      save 
      DIMENSION P00(NS),P10(NS),P11(NS),P20(NS),P21(NS),Q21(NS),
     *          COA(IM,JM,NT,NS),COP(IM,JM,NT,NS),FRQ(NT,NS),
     *          A(NT,NT),IND(NT),XB(NT),YB(NT),HBR(NT,NS),
     *          C(NT,NS),S(NT,NS),
     *          X(NT),Y(NT)
C     DIMENSION U(IM,JM,NT,NS),V(IM,JM,NT,NS)
      DIMENSION INDX(NC,NS,4),AMP(NC,NS),FREQ(NC,NS),
     *          PHA(NC,NS),CT(NC,NS),ST(NC,NS),DPD(4)
      REAL TIME,h1(IM,JM,NT,NS),H2(IM,JM,NT,NS),
     *          ph1(nt,ns),ph2(nt,ns),H(IM,JM),u,v
c 
      DATA P00/ 0.0298D0, 0.0200D0/, P10/0.1408D0,0.0905D0/,
     *     P11/-0.0805D0,-0.0638D0/, P20/0.6002D0,0.3476D0/,
     *     P21/-0.3025D0,-0.1645D0/, Q21/0.1517D0,0.0923D0/
      DATA DPD/13.1763965D0,0.9856473D0,0.1114041D0,0.0529539D0/
c
      COMMON/OTD/u(im,jm,nt,ns),v(im,jm,nt,ns)  
	logical readneed
        data readneed/.true./


	if(readneed) then
          write(*,*) "Read Topo" 
	  open(61,file='gtdtdu40e',form='formatted')
	  read(61,120)h
120	format(8(1pe10.3))
	close(61)
	readneed=.FALSE.
      DT=2.D0
      RAD = DATAN(1.D0)/45.D0
      TWOPI = DATAN(1.D0)*8.D0
      PI=TWOPI*0.5D0
      DPLD  = 36.D1 - DPD(1) + DPD(2)

      OPEN (8,file='doodson.dat',form='formatted')
      FDPD  = 0.D0
C     WRITE(6,150)
C 150 FORMAT(/'   Semidiurnal & diurnal tidal potential terms:'/)
      DO 180 M=1,2
      FDPD = FDPD + DPLD
      DO 170 L=1,30
      READ (8,*) (INDX(L,M,N),N=1,4), AMP(L,M), PHA(L,M)
      THETA = FDPD
      DO 160 N=1,4
      THETA = THETA + REAL(INDX(L,M,N))*DPD(N)
  160 CONTINUE
      FREQ(L,M) = THETA
      OMEGAT  = 2.D0*THETA*RAD
      CT(L,M) = 2.D0*DCOS(OMEGAT)
      ST(L,M) = 2.D0*DSIN(OMEGAT)
C     WRITE(6,502) L,M,(INDX(L,M,N),N=1,4),FREQ(L,M),AMP(L,M),
C    *             PHA(L,M)
C 502 FORMAT(1X,I3,I4,2X,4I2,F12.6,2(F10.5,F10.1,5X))
  170 CONTINUE
      WRITE(6,*) ' '
  180 CONTINUE

      FRQ(1,1)=FREQ(20,1)
      FRQ(2,1)=FREQ( 9,1)
      FRQ(3,1)=FREQ(18,1)
      FRQ(1,2)=FREQ(15,2)
      FRQ(2,2)=FREQ(22,2)
      FRQ(3,2)=FREQ( 8,2)

      HBR(1,1)=AMP(20,1)
      HBR(2,1)=AMP( 9,1)
      HBR(3,1)=AMP(18,1)
      HBR(1,2)=AMP(15,2)
      HBR(2,2)=AMP(22,2)
      HBR(3,2)=AMP( 8,2)

      C(1,1)=CT(20,1)
      C(2,1)=CT( 9,1)
      C(3,1)=CT(18,1)
      C(1,2)=CT(15,2)
      C(2,2)=CT(22,2)
      C(3,2)=CT( 8,2)
      S(1,1)=ST(20,1)
      S(2,1)=ST( 9,1)
      S(3,1)=ST(18,1)
      S(1,2)=ST(15,2)
      S(2,2)=ST(22,2)
      S(3,2)=ST( 8,2)
C END if NEEDREAD 
C     ENDIF



C     OPEN(60,FILE='/archive/students/liaom/gtd.ortd.all',status='old')
C     READ (60,120) U,V
C     CLOSE(60)
      write(*,*)'Processed ortd file!'

C     CONVERSION OF U,V >> COA,COP

      DO 250 M=1,NS
      DO 250 L=1,NT
      DO 240 J=1,JM
      DO 240 I=1,IM
      COA(I,J,L,M)= 0.
      COP(I,J,L,M)= 0.
      IF (H(I,J).LT.2.) GO TO 240
      XA=(P00(M))*U(I,J,1,M)+
     1   (P10(M)+P11(M)*C(L,M))*U(I,J,2,M)+
     2   (P20(M)+P21(M)*C(L,M)+Q21(M)*S(L,M))*U(I,J,3,M)
      YA=(P00(M))*V(I,J,1,M)+
     1   (P10(M)+P11(M)*C(L,M))*V(I,J,2,M)+
     2   (P20(M)+P21(M)*C(L,M)+Q21(M)*S(L,M))*V(I,J,3,M)
      IF((ABS(XA).LT.SMALL).AND.(XA.GE.0.0))  XA = SMALL 
      IF((ABS(XA).LT.SMALL).AND.(XA.LT.0.0))  XA = -1.0D0*SMALL 
      COA(I,J,L,M)=DSQRT(XA**2+YA**2)*HBR(L,M)
      COP(I,J,L,M)=-DATAN2(YA,XA)/RAD
      IF(M.EQ.1) COP(I,J,L,M)=COP(I,J,L,M)+PI/RAD
      IF(COP(I,J,L,M).LT.0.) COP(I,J,L,M)=COP(I,J,L,M)+2.*PI/RAD
  240 CONTINUE
  250 CONTINue
      
      DO 230 M=1,NS
      DO 230 L=1,NT
      DO 220 J=1,JM
      DO 220 I=1,IM
        h1(i,j,l,m)=coa(i,j,l,m)*cos(cop(i,j,l,m)*pi/180.)
        h2(i,j,l,m)=coa(i,j,l,m)*sin(cop(i,j,l,m)*pi/180.)
  220 CONTINUE
  230 CONTINUE
  
C END if NEEDREAD 
      END IF 
      call cotid(REAL(alon),REAL(alat),im,jm,nt,h1,h2,ns,ph1,ph2)

      
1000  RETURN    
      END

      SUBROUTINE CRAMER(A,NT,B,X)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(NT,NT),B(NT),X(NT)

C     VALID ONLY FOR NT=3

      DENM=A(1,1)*(A(2,2)*A(3,3)-A(3,2)*A(2,3))
     1    -A(1,2)*(A(2,1)*A(3,3)-A(3,1)*A(2,3))
     2    +A(1,3)*(A(2,1)*A(3,2)-A(3,1)*A(2,2))
      X(1)=B(1  )*(A(2,2)*A(3,3)-A(3,2)*A(2,3))
     1    -A(1,2)*(B(2  )*A(3,3)-B(3  )*A(2,3))
     2    +A(1,3)*(B(2  )*A(3,2)-B(3  )*A(2,2))
      X(2)=A(1,1)*(B(2  )*A(3,3)-B(3  )*A(2,3))
     1    -B(1  )*(A(2,1)*A(3,3)-A(3,1)*A(2,3))
     2    +A(1,3)*(A(2,1)*B(3  )-A(3,1)*B(2  ))
      X(3)=A(1,1)*(A(2,2)*B(3  )-A(3,2)*B(2  ))
     1    -A(1,2)*(A(2,1)*B(3  )-A(3,1)*B(2  ))
     2    +B(1  )*(A(2,1)*A(3,2)-A(3,1)*A(2,2))
      X(1)=X(1)/DENM
      X(2)=X(2)/DENM
      X(3)=X(3)/DENM

      RETURN
      END
      subroutine ludcmp(a,n,np,indx,d)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      parameter (nmax=100,tiny=1.0e-20)
      dimension a(np,np),indx(n),vv(nmax)
      d=1.
      do 12 i=1,n
        aamax=0.
        do 11 j=1,n
          if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
11      continue
        if (aamax.eq.0.) pause 'singular matrix.'
        vv(i)=1./aamax
12    continue
      do 19 j=1,n
        if (j.gt.1) then
          do 14 i=1,j-1
            sum=a(i,j)
            if (i.gt.1)then
              do 13 k=1,i-1
                sum=sum-a(i,k)*a(k,j)
13            continue
              a(i,j)=sum
            endif
14        continue
        endif
        aamax=0.
        do 16 i=j,n
          sum=a(i,j)
          if (j.gt.1)then
            do 15 k=1,j-1
              sum=sum-a(i,k)*a(k,j)
15          continue
            a(i,j)=sum
          endif
          dum=vv(i)*abs(sum)
          if (dum.ge.aamax) then
            imax=i
            aamax=dum
          endif
16      continue
        if (j.ne.imax)then
          do 17 k=1,n
            dum=a(imax,k)
            a(imax,k)=a(j,k)
            a(j,k)=dum
17        continue
          d=-d
          vv(imax)=vv(j)
        endif
        indx(j)=imax
        if(j.ne.n)then
          if(a(j,j).eq.0.)a(j,j)=tiny
          dum=1./a(j,j)
          do 18 i=j+1,n
            a(i,j)=a(i,j)*dum
18        continue
        endif
19    continue
      if(a(n,n).eq.0.)a(n,n)=tiny
      return
      end

      subroutine lubksb(a,n,np,indx,b)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      dimension a(np,np),indx(n),b(n)
      ii=0
      do 12 i=1,n
        ll=indx(i)
        sum=b(ll)
        b(ll)=b(i)
        if (ii.ne.0)then
          do 11 j=ii,i-1
            sum=sum-a(i,j)*b(j)
11        continue
        else if (sum.ne.0.) then
          ii=i
        endif
        b(i)=sum
12    continue
      do 14 i=n,1,-1
        sum=b(i)
        if(i.lt.n)then
          do 13 j=i+1,n
            sum=sum-a(i,j)*b(j)
13        continue
        endif
        b(i)=sum/a(i,i)
14    continue
      return
      end


