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	
