C ********************************************
C * Tide Calculations *
C * *
C * Program Copyright (C) 1996 Alan Bain *
C ********************************************
C Disclaimer:
C This program was written for my own interest, to learn
C how tidal predictions can be made. I don't think there
C are any errors, but I take no responsibility for any
C damage, direct or indirect caused by the use of results
C from this program
C
C This program may be freely distributed, however I would
C appreciate to know if you find it useful,
C
C alanb@chiark.greenend.org.uk
C Instructions
C The program is written in Fortran 77 and should compile on
C most systems. Please let me know at the above address if
C you have problems compiling it.
C To make predictions you need this program and a data file
C for the port of interest. I include a Liverpool one using
C data obtained from the Admiralty listing of tidal harmonics
C The constituent breakdown comes from a paper of Mohamed Amin.
C
C The port data should be kept in the file port.data.
C To compile on a standard unix system the following
C should worka (assuming you save this file as tide.f)
C
C f77 -o tide tide.f
C ./tide
C
C Improvements
C ------------
C I should like to make the following improvements, when
C I have time :-
C (a) Use of standard stationary point technique to
C find high and low waters on a given day.
C (this needs dh/dt=0 as condition)
C (b) Improved output handling
C (c) Tidal differences support.
C (d) Add leap centuries support.
C (e) Update the ephemeris -- Brown is rather old.
C Please mail me if you wish to suggest any other improvements.
C
C Warning -- the ephemeris is only accurate for ports on W coast
C of England -- I shall improve this in due course.
C
C Format of port data file
C (No. consituents 100 max.)
C (Chart Datum)
C N times (r,a,b,c,d,e,phi,sig,h,g)
C note -- r,a,b,c,d,e,phi,sig are constant for each constituent,
C only h and g need to be set for the individual port.
C see Doodson, Cartwright or Amin for tables of r,a,b,c,d,e,phi
C for several distinct decompositions. The _first_ constituent
C should (following common practice) be M2.
PROGRAM CALC_TIDE
IMPLICIT NONE
INTEGER NH
PARAMETER (NH=100)
REAL SIG(NH),G(NH),Hz(NH),R(NH),A(NH),B(NH),PHI(NH)
REAL C(NH),D(NH),E(NH)
REAL s,h,p,n,pp,pi,dtr
REAL V,fM2,uM2,T,tau,Long
REAL CDATUM
INTEGER nday,ndays,year,i,j,nmax
REAL time,height,contrib
integer month(12),nmonth
DATA Pi/3.14159245368979/
C Month Lengths
MONTH(1)=31
MONTH(2)=28
MONTH(3)=31
MONTH(4)=30
MONTH(5)=31
MONTH(6)=30
MONTH(7)=31
MONTH(8)=31
MONTH(9)=30
MONTH(10)=31
MONTH(11)=30
MONTH(12)=31
DTR=Pi/180.0
C This is the data file
WRITE(*,*)'Tide Calculator -- Alan Bain'
WRITE(*,*)'Any comments, please mail alanb@chiark.greenend.org.uk'
OPEN(20,'port.data',status='old',form='formatted')
C Mean solar time is (time) and starts at midnight GMT
READ(20,*)NMAX
READ(20,*)CDATUM
height=0.0
write(*,*)'Year?'
read(*,*)year
write(*,*)'Month?'
read(*,*)nmonth
write(*,*)'Day ?'
read(*,*)nday
write(*,*)'Time (hours from Midnight)?'
read(*,*)time
write(*,*)'Calculating Tides for Port'
write(*,*)'Year : ',year
write(*,*)'Day : ',nday
write(*,*)'Month : ',nmonth
write(*,*)'Time : ',time
Do I=1,NMAX
READ(20,*)R(I),A(i),B(i),C(i),D(i),E(i),PHI(I),
& SIG(I),Hz(I),G(I)
C This routine needs year and nday (days into the year).
ndays=0
Do J=1,(YEAR-1900)
IF ((J+1900)/4*4.eq.(j+1900)) then
ndays=ndays+366
month(2)=29
else
ndays=ndays+365
month(2)=28
endif
enddo
do j=1,(nmonth-1)
ndays=ndays+month(j)
enddo
C No leap centuries are allowed
ndays=ndays+(nday-1)
T=ndays/36525.0
write(*,*)T
C Convert mean solar days to Julian Centuries
C From Doodson 1921a
C Based on the ephemeris of Brown 1907
s=277.0248+481267.8906*T+0.0020*T*T
h=280.1895+36000.7689*T+0.0003*T*T
p=334.3853+4069.0340*T-0.0103*T*T
N=100.8432+1934.1420*T-0.0021*T*T
pp=281.2209+1.7192*T+0.0005*T*T
C Here T is in julian centuries since Midnight
C at Grenwich on the 1st Jan 1900
C Note a julian century is 26,525 mean solar days
C Now we know the astronomy at the required time
C V=r(i)*tau+a(i)*s+b(i)*h+c(i)*p+d(i)*N+e(i)*pp+phi(i)
V=(a(i)-r(i))*s+(b(i)+r(i))*h+c(i)*p+d(i)*N+e(i)*pp+
& phi(i)
C corrections to M2
uM2=-2.1*sin(N*DTR)
fM2=1.0-0.037*cos(N*DTR)
C Now to find the height of water contributed by this
C component calculate H cos (V+sig(i)*t+u-G)
if (i.ne.1) then
uM2=0.0
fM2=1.0
endif
contrib=Hz(i)*fM2*cos((V+sig(i)*time+uM2-G(i)+phi(i))*DTR)
height=height+contrib
enddo
height=height+cdatum
C Datum for the chart (usually obtained from the Admiralty Handbook)
write(*,*)'Height of Water (m) ',height
end