      PROGRAM TIDE

C  THIS PROGRAM MUST BE RUN WITH DOUBLE-PRECISION REAL NUMBERS

      DIMENSION PHIX(0:1000),PHIY(0:1000)

C  MASS OF EARTH, MOON AND EARTH-MOON DISTANCE
      DATA EM1,EM2,R12/5.976E24,7.35E22,3.84E8/

C  RADIUS OF EARTH AND MEAN DEPTH OF OCEAN
      DATA RE,D/6.378E6,2132/  

C  FIND POTENTIAL ALONG X AXIS WITHIN REGION AROUND OCEAN SURFACE

      DO 1 I=0,1000
        XI=0.01*I-5.00
        PHIX(I)=-EM1/(RE+XI)-EM2/(R12-RE-XI)+EM2*(RE+XI)/R12**2
        PHIY(I)=-EM1/(RE+XI)-EM2/R12
    1 CONTINUE

C  FIND MINIMUM DIFFERENCE OF  POTENTIALS AT HEIGHT H ABOVE SURFACE IN
C  MOON DIRECTION AND DEPTH 0.5*H PERPENDICULAR TO MOON DIRECTION.
C  THIS WILL GIVE CORRECT VOLUME AND INDICATE CORRECT MAXIMUM TIDE HEIGHT.

      DIFF=1E20
      DO 3 I=0,250
        ERR=ABS(PHIX(500+2*I)-PHIY(500-I))
        IF(ERR.LT.DIFF)THEN
                DIFF=ERR
                HEIGHT=0.02*I
        ENDIF
    3 CONTINUE
      WRITE(6,100)HEIGHT
  100 FORMAT(15H HEIGHT OF TIDE, F6.2,7H METRES)

C  FOR GRAPHICAL OUTPUT CREATE FILE GIVING CORRESPONDING POTENTIALS
C  AT HEIGHT H AT POINT N AND DEPTH 0.5*H AT POINT S.
C  POTENTIALS ARE EXPRESSED IN UNITS OF 10**-10 TIMES THE DIFFERENCE
C  THE POTENTIAL AND THAT AT N FOR THE UNDISTURBED SURFACE.

      OPEN(UNIT=10,FILE='TIDE.DAT')
      DO 4 I=0,250
        HEIGHT=0.02*I
        F1=(PHIX(500+2*I)-PHIX(500))/1E10
        F2=(PHIY(500-I)-PHIX(500))/1E10
        WRITE(10,*)HEIGHT,F1,F2
    4 CONTINUE

      STOP
      END
      
